Distributed-Run-Time-Resource-Manager/RCCE_V2.0/apps/NPB/LU/jacu.c

// 
// Copyright 2010 Intel Corporation
// 
//    Licensed under the Apache License, Version 2.0 (the "License");
//    you may not use this file except in compliance with the License.
//    You may obtain a copy of the License at
// 
//        http://www.apache.org/licenses/LICENSE-2.0
// 
//    Unless required by applicable law or agreed to in writing, software
//    distributed under the License is distributed on an "AS IS" BASIS,
//    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//    See the License for the specific language governing permissions and
//    limitations under the License.
// 
#include "applu_share.h"
#include "applu_macros.h"

void jacu(int k) {

//   compute the upper triangular part of the jacobian matrix

//  local variables

  int i, j;
  double  r43;
  double  c1345;
  double  c34;
  double  tmp1, tmp2, tmp3;

  r43 = ( 4.0 / 3.0 );
  c1345 = c1 * c3 * c4 * c5;
  c34 = c3 * c4;

  for (j = jst; j<= jend; j++) 
    for (i = ist; i<=iend; i++) {

//c---------------------------------------------------------------------
//c   form the block diagonal
//c---------------------------------------------------------------------
               tmp1 = 1.0 / u(1,i,j,k);
               tmp2 = tmp1 * tmp1;
               tmp3 = tmp1 * tmp2;

               tmp1 = 1.0 / u(1,i,j,k);
               tmp2 = tmp1 * tmp1;
               tmp3 = tmp1 * tmp2;

               d(1,1,i,j) =  1.0
                             + dt * 2.0 * (   tx1 * dx1
                                                + ty1 * dy1
                                                + tz1 * dz1 );
               d(1,2,i,j) =  0.0;
               d(1,3,i,j) =  0.0;
               d(1,4,i,j) =  0.0;
               d(1,5,i,j) =  0.0;

               d(2,1,i,j) =  dt * 2.0
                * (  tx1 * ( - r43 * c34 * tmp2 * u(2,i,j,k) )
                   + ty1 * ( -       c34 * tmp2 * u(2,i,j,k) )
                   + tz1 * ( -       c34 * tmp2 * u(2,i,j,k) ) );
               d(2,2,i,j) =  1.0
                + dt * 2.0 
                * (  tx1 * r43 * c34 * tmp1
                   + ty1 *       c34 * tmp1
                   + tz1 *       c34 * tmp1 )
                + dt * 2.0 * (   tx1 * dx2
                                   + ty1 * dy2
                                   + tz1 * dz2  );
               d(2,3,i,j) = 0.0;
               d(2,4,i,j) = 0.0;
               d(2,5,i,j) = 0.0;

               d(3,1,i,j) = dt * 2.0
            * (  tx1 * ( -       c34 * tmp2 * u(3,i,j,k) )
               + ty1 * ( - r43 * c34 * tmp2 * u(3,i,j,k) )
               + tz1 * ( -       c34 * tmp2 * u(3,i,j,k) ) );
               d(3,2,i,j) = 0.0;
               d(3,3,i,j) = 1.0
               + dt * 2.0
                    * (  tx1 *       c34 * tmp1
                       + ty1 * r43 * c34 * tmp1
                       + tz1 *       c34 * tmp1 )
               + dt * 2.0 * (  tx1 * dx3
                                 + ty1 * dy3
                                 + tz1 * dz3 );
               d(3,4,i,j) = 0.0;
               d(3,5,i,j) = 0.0;

               d(4,1,i,j) = dt * 2.0
            * (  tx1 * ( -       c34 * tmp2 * u(4,i,j,k) )
               + ty1 * ( -       c34 * tmp2 * u(4,i,j,k) )
               + tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k) ) );
               d(4,2,i,j) = 0.0;
               d(4,3,i,j) = 0.0;
               d(4,4,i,j) = 1.0
               + dt * 2.0
                    * (  tx1 *       c34 * tmp1
                       + ty1 *       c34 * tmp1
                       + tz1 * r43 * c34 * tmp1 )
               + dt * 2.0 * (  tx1 * dx4
                                 + ty1 * dy4
                                 + tz1 * dz4 );
               d(4,5,i,j) = 0.0;

               d(5,1,i,j) = dt * 2.0
       * ( tx1 * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i,j,k) *u(2,i,j,k) )
                   - ( c34 - c1345 ) * tmp3 * ( u(3,i,j,k) * u(3,i,j,k) )
                   - ( c34 - c1345 ) * tmp3 * ( u(4,i,j,k) * u(4,i,j,k) )
                   - ( c1345 ) * tmp2 * u(5,i,j,k) )
         + ty1 * ( - ( c34 - c1345 ) * tmp3 * ( u(2,i,j,k) * u(2,i,j,k) )
                   - ( r43*c34 - c1345 ) * tmp3 * ( u(3,i,j,k) * u(3,i,j,k) )
                   - ( c34 - c1345 ) * tmp3 * ( u(4,i,j,k) * u(4,i,j,k) )
                   - ( c1345 ) * tmp2 * u(5,i,j,k) )
         + tz1 * ( - ( c34 - c1345 ) * tmp3 * ( u(2,i,j,k) * u(2,i,j,k) )
                   - ( c34 - c1345 ) * tmp3 * ( u(3,i,j,k) * u(3,i,j,k) )
                   - ( r43*c34 - c1345 ) * tmp3 * ( u(4,i,j,k) *u(4,i,j,k) )
                   - ( c1345 ) * tmp2 * u(5,i,j,k) ) );
               d(5,2,i,j) = dt * 2.0
       * ( tx1 * ( r43*c34 - c1345 ) * tmp2 * u(2,i,j,k)
         + ty1 * (     c34 - c1345 ) * tmp2 * u(2,i,j,k)
         + tz1 * (     c34 - c1345 ) * tmp2 * u(2,i,j,k) );
               d(5,3,i,j) = dt * 2.0
       * ( tx1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k)
         + ty1 * ( r43*c34 -c1345 ) * tmp2 * u(3,i,j,k)
         + tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k) );
               d(5,4,i,j) = dt * 2.0
       * ( tx1 * ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
         + ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
         + tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k) );
               d(5,5,i,j) = 1.0
         + dt * 2.0 * ( tx1 * c1345 * tmp1
                          + ty1 * c1345 * tmp1
                          + tz1 * c1345 * tmp1 )
         + dt * 2.0 * (  tx1 * dx5
                          +  ty1 * dy5
                          +  tz1 * dz5 );

//c---------------------------------------------------------------------
//c   form the first block sub-diagonal
//c---------------------------------------------------------------------
               tmp1 = 1.0 / u(1,i+1,j,k);
               tmp2 = tmp1 * tmp1;
               tmp3 = tmp1 * tmp2;

               a(1,1,i,j) = - dt * tx1 * dx1;
               a(1,2,i,j) =   dt * tx2;
               a(1,3,i,j) =   0.0;
               a(1,4,i,j) =   0.0;
               a(1,5,i,j) =   0.0;

               a(2,1,i,j) =  dt * tx2
                * ( - ( u(2,i+1,j,k) * tmp1 ) *( u(2,i+1,j,k) * tmp1 )
           + c2 * 0.50 * (  u(2,i+1,j,k) * u(2,i+1,j,k)
                              + u(3,i+1,j,k) * u(3,i+1,j,k)
                              + u(4,i+1,j,k) * u(4,i+1,j,k) ) * tmp2 )
                - dt * tx1 * ( - r43 * c34 * tmp2 * u(2,i+1,j,k) );
               a(2,2,i,j) =  dt * tx2
                * ( ( 2.0 - c2 ) * ( u(2,i+1,j,k) * tmp1 ) )
                - dt * tx1 * ( r43 * c34 * tmp1 )
                - dt * tx1 * dx2;
               a(2,3,i,j) =  dt * tx2
                    * ( - c2 * ( u(3,i+1,j,k) * tmp1 ) );
               a(2,4,i,j) =  dt * tx2
                    * ( - c2 * ( u(4,i+1,j,k) * tmp1 ) );
               a(2,5,i,j) =  dt * tx2 * c2 ;

               a(3,1,i,j) =  dt * tx2
                    * ( - ( u(2,i+1,j,k) * u(3,i+1,j,k) ) * tmp2 )
               - dt * tx1 * ( - c34 * tmp2 * u(3,i+1,j,k) );
               a(3,2,i,j) =  dt * tx2 * ( u(3,i+1,j,k) * tmp1 );
               a(3,3,i,j) =  dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
                - dt * tx1 * ( c34 * tmp1 )
                - dt * tx1 * dx3;
               a(3,4,i,j) = 0.0;
               a(3,5,i,j) = 0.0;

               a(4,1,i,j) = dt * tx2
                * ( - ( u(2,i+1,j,k)*u(4,i+1,j,k) ) * tmp2 )
                - dt * tx1 * ( - c34 * tmp2 * u(4,i+1,j,k) );
               a(4,2,i,j) = dt * tx2 * ( u(4,i+1,j,k) * tmp1 );
               a(4,3,i,j) = 0.0;
               a(4,4,i,j) = dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
                - dt * tx1 * ( c34 * tmp1 )
                - dt * tx1 * dx4;
               a(4,5,i,j) = 0.0;

               a(5,1,i,j) = dt * tx2
                * ( ( c2 * (  u(2,i+1,j,k) * u(2,i+1,j,k)
                            + u(3,i+1,j,k) * u(3,i+1,j,k)
                            + u(4,i+1,j,k) * u(4,i+1,j,k) ) * tmp2
                    - c1 * ( u(5,i+1,j,k) * tmp1 ) )
                * ( u(2,i+1,j,k) * tmp1 ) )
                - dt * tx1
                * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i+1,j,k)*u(2,i+1,j,k) )
                    - (     c34 - c1345 ) * tmp3 * ( u(3,i+1,j,k)*u(3,i+1,j,k) )
                    - (     c34 - c1345 ) * tmp3 * ( u(4,i+1,j,k)*u(4,i+1,j,k) )
                    - c1345 * tmp2 * u(5,i+1,j,k) );
               a(5,2,i,j) = dt * tx2
                * ( c1 * ( u(5,i+1,j,k) * tmp1 )
                   - 0.50 * c2
                   * ( (  3.0*u(2,i+1,j,k)*u(2,i+1,j,k)
                        + u(3,i+1,j,k)*u(3,i+1,j,k)
                        + u(4,i+1,j,k)*u(4,i+1,j,k) ) * tmp2 ) )
                 - dt * tx1
                 * ( r43*c34 - c1345 ) * tmp2 * u(2,i+1,j,k);
               a(5,3,i,j) = dt * tx2
                 * ( - c2 * ( u(3,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
                 - dt * tx1
                 * (  c34 - c1345 ) * tmp2 * u(3,i+1,j,k);
               a(5,4,i,j) = dt * tx2
                 * ( - c2 * ( u(4,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
                 - dt * tx1
                 * (  c34 - c1345 ) * tmp2 * u(4,i+1,j,k);
               a(5,5,i,j) = dt * tx2
                 * ( c1 * ( u(2,i+1,j,k) * tmp1 ) )
                 - dt * tx1 * c1345 * tmp1
                 - dt * tx1 * dx5;

//c---------------------------------------------------------------------
//c   form the second block sub-diagonal
//c---------------------------------------------------------------------
               tmp1 = 1.0 / u(1,i,j+1,k);
               tmp2 = tmp1 * tmp1;
               tmp3 = tmp1 * tmp2;

               b(1,1,i,j) = - dt * ty1 * dy1;
               b(1,2,i,j) =   0.0;
               b(1,3,i,j) =  dt * ty2;
               b(1,4,i,j) =   0.0;
               b(1,5,i,j) =   0.0;

               b(2,1,i,j) =  dt * ty2
                 * ( - ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
                 - dt * ty1 * ( - c34 * tmp2 * u(2,i,j+1,k) );
               b(2,2,i,j) =  dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
                - dt * ty1 * ( c34 * tmp1 )
                - dt * ty1 * dy2;
               b(2,3,i,j) =  dt * ty2 * ( u(2,i,j+1,k) * tmp1 );
               b(2,4,i,j) = 0.0;
               b(2,5,i,j) = 0.0;

               b(3,1,i,j) =  dt * ty2
                 * ( - ( u(3,i,j+1,k) * tmp1 ) * ( u(3,i,j+1,k) * tmp1 )
            + 0.50 * c2 * ( (  u(2,i,j+1,k) * u(2,i,j+1,k)
                                 + u(3,i,j+1,k) * u(3,i,j+1,k)
                                 + u(4,i,j+1,k) * u(4,i,j+1,k) )
                                * tmp2 ) )
             - dt * ty1 * ( - r43 * c34 * tmp2 * u(3,i,j+1,k) );
               b(3,2,i,j) =  dt * ty2
                         * ( - c2 * ( u(2,i,j+1,k) * tmp1 ) );
               b(3,3,i,j) =  dt * ty2 * ( ( 2.0 - c2 )
                         * ( u(3,i,j+1,k) * tmp1 ) )
             - dt * ty1 * ( r43 * c34 * tmp1 )
             - dt * ty1 * dy3;
               b(3,4,i,j) =  dt * ty2
                         * ( - c2 * ( u(4,i,j+1,k) * tmp1 ) );
               b(3,5,i,j) =  dt * ty2 * c2;

               b(4,1,i,j) =  dt * ty2
                    * ( - ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
             - dt * ty1 * ( - c34 * tmp2 * u(4,i,j+1,k) );
               b(4,2,i,j) = 0.0;
               b(4,3,i,j) =  dt * ty2 * ( u(4,i,j+1,k) * tmp1 );
               b(4,4,i,j) =  dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
                              - dt * ty1 * ( c34 * tmp1 )
                              - dt * ty1 * dy4;
               b(4,5,i,j) = 0.0;

               b(5,1,i,j) =  dt * ty2
                * ( ( c2 * (  u(2,i,j+1,k) * u(2,i,j+1,k)
                            + u(3,i,j+1,k) * u(3,i,j+1,k)
                            + u(4,i,j+1,k) * u(4,i,j+1,k) ) * tmp2
                     - c1 * ( u(5,i,j+1,k) * tmp1 ) )
                * ( u(3,i,j+1,k) * tmp1 ) )
                - dt * ty1
                * ( - (     c34 - c1345 )*tmp3*(u(2,i,j+1,k)*u(2,i,j+1,k))
                    - ( r43*c34 - c1345 )*tmp3*(u(3,i,j+1,k)*u(3,i,j+1,k))
                    - (     c34 - c1345 )*tmp3*(u(4,i,j+1,k)*u(4,i,j+1,k))
                    - c1345*tmp2*u(5,i,j+1,k) );
               b(5,2,i,j) =  dt * ty2
                * ( - c2 * ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
                - dt * ty1
                * ( c34 - c1345 ) * tmp2 * u(2,i,j+1,k);
               b(5,3,i,j) =  dt * ty2
                * ( c1 * ( u(5,i,j+1,k) * tmp1 )
                - 0.50 * c2 
                * ( (  u(2,i,j+1,k)*u(2,i,j+1,k)
                     + 3.0 * u(3,i,j+1,k)*u(3,i,j+1,k)
                     + u(4,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 ) )
                - dt * ty1
                * ( r43*c34 - c1345 ) * tmp2 * u(3,i,j+1,k);
               b(5,4,i,j) =  dt * ty2
                * ( - c2 * ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
                - dt * ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j+1,k);
               b(5,5,i,j) =  dt * ty2
                * ( c1 * ( u(3,i,j+1,k) * tmp1 ) )
                - dt * ty1 * c1345 * tmp1
                - dt * ty1 * dy5;

//c---------------------------------------------------------------------
//c   form the third block sub-diagonal
//c---------------------------------------------------------------------
               tmp1 = 1.0 / u(1,i,j,k+1);
               tmp2 = tmp1 * tmp1;
               tmp3 = tmp1 * tmp2;

               c(1,1,i,j) = - dt * tz1 * dz1;
               c(1,2,i,j) =   0.0;
               c(1,3,i,j) =   0.0;
               c(1,4,i,j) = dt * tz2;
               c(1,5,i,j) =   0.0;

               c(2,1,i,j) = dt * tz2
                 * ( - ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
                 - dt * tz1 * ( - c34 * tmp2 * u(2,i,j,k+1) );
               c(2,2,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
                 - dt * tz1 * c34 * tmp1
                 - dt * tz1 * dz2 ;
               c(2,3,i,j) = 0.0;
               c(2,4,i,j) = dt * tz2 * ( u(2,i,j,k+1) * tmp1 );
               c(2,5,i,j) = 0.0;

               c(3,1,i,j) = dt * tz2
                 * ( - ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
                 - dt * tz1 * ( - c34 * tmp2 * u(3,i,j,k+1) );
               c(3,2,i,j) = 0.0;
               c(3,3,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
                 - dt * tz1 * ( c34 * tmp1 )
                 - dt * tz1 * dz3;
               c(3,4,i,j) = dt * tz2 * ( u(3,i,j,k+1) * tmp1 );
               c(3,5,i,j) = 0.0;

               c(4,1,i,j) = dt * tz2
              * ( - ( u(4,i,j,k+1) * tmp1 ) *( u(4,i,j,k+1) * tmp1 )
                  + 0.50 * c2
                  * ( ( u(2,i,j,k+1) * u(2,i,j,k+1)
                      + u(3,i,j,k+1) * u(3,i,j,k+1)
                      + u(4,i,j,k+1) * u(4,i,j,k+1) ) * tmp2 ) )
              - dt * tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k+1) );
               c(4,2,i,j) = dt * tz2
                   * ( - c2 * ( u(2,i,j,k+1) * tmp1 ) );
               c(4,3,i,j) = dt * tz2
                   * ( - c2 * ( u(3,i,j,k+1) * tmp1 ) );
               c(4,4,i,j) = dt * tz2 * ( 2.0 - c2 )
                   * ( u(4,i,j,k+1) * tmp1 )
                   - dt * tz1 * ( r43 * c34 * tmp1 )
                   - dt * tz1 * dz4;
               c(4,5,i,j) = dt * tz2 * c2;

               c(5,1,i,j) = dt * tz2
           * ( ( c2 * (  u(2,i,j,k+1) * u(2,i,j,k+1)
                       + u(3,i,j,k+1) * u(3,i,j,k+1)
                       + u(4,i,j,k+1) * u(4,i,j,k+1) ) * tmp2
             - c1 * ( u(5,i,j,k+1) * tmp1 ) )
                  * ( u(4,i,j,k+1) * tmp1 ) )
             - dt * tz1
             * ( - ( c34 - c1345 ) * tmp3 * (u(2,i,j,k+1)*u(2,i,j,k+1))
                 - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k+1)*u(3,i,j,k+1))
                 - ( r43*c34 - c1345 )* tmp3 * (u(4,i,j,k+1)*u(4,i,j,k+1))
                - c1345 * tmp2 * u(5,i,j,k+1) );
               c(5,2,i,j) = dt * tz2
             * ( - c2 * ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
             - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k+1);
               c(5,3,i,j) = dt * tz2
             * ( - c2 * ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
             - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k+1);
               c(5,4,i,j) = dt * tz2
             * ( c1 * ( u(5,i,j,k+1) * tmp1 )
             - 0.50 * c2
             * ( (  u(2,i,j,k+1)*u(2,i,j,k+1)
                  + u(3,i,j,k+1)*u(3,i,j,k+1)
                  + 3.0*u(4,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 ) )
             - dt * tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k+1);
               c(5,5,i,j) = dt * tz2
             * ( c1 * ( u(4,i,j,k+1) * tmp1 ) )
             - dt * tz1 * c1345 * tmp1
             - dt * tz1 * dz5;

      }
      return;
}