MOAB: Mesh Oriented datABase  (version 5.4.0)
tensor.c
Go to the documentation of this file.
00001 #include "moab/FindPtFuncs.h"
00002
00003 /*--------------------------------------------------------------------------
00004    Matrix-Matrix Multiplication
00005
00006    mxm_ab (A,na,B,nb,C,nc) :
00007       gives C = A B where A is na x nb, B is nb x nc, C is na x nc
00008       a := r | c   to indicate A is in row- or column- major format
00009       b := r | c   to indicate B is in row- or column- major format
00010       C is always column-major
00011   --------------------------------------------------------------------------*/
00012
00013 static void mxm_cc( const realType* A, unsigned na, const realType* B, unsigned nb, realType* C, unsigned nc )
00014 {
00015     unsigned i, j, k;
00016     realType* Ccol       = C;
00017     const realType* Bcol = B;
00018     for( j = 0; j < nc; ++j, Ccol += na, Bcol += nb )
00019     {
00020         const realType* Acol = A;
00021         for( i = 0; i < na; ++i )
00022             Ccol[i] = 0;
00023         for( k = 0; k < nb; ++k, Acol += na )
00024             for( i = 0; i < na; ++i )
00025                 Ccol[i] += Acol[i] * Bcol[k];
00026     }
00027 }
00028
00029 static void mxm_rc( const realType* A, unsigned na, const realType* B, unsigned nb, realType* C, unsigned nc )
00030 {
00031     unsigned i, j, k;
00032     realType* Ccol       = C;
00033     const realType* Bcol = B;
00034     for( j = 0; j < nc; ++j, Ccol += na, Bcol += nb )
00035     {
00036         const realType* Arow = A;
00037         for( i = 0; i < na; ++i, Arow += nb )
00038         {
00039             Ccol[i] = 0;
00040             for( k = 0; k < nb; ++k )
00041                 Ccol[i] += Arow[k] * Bcol[k];
00042         }
00043     }
00044 }
00045
00046 static void mxm_cr( const realType* A, unsigned na, const realType* B, unsigned nb, realType* C, unsigned nc )
00047 {
00048     unsigned i, j, k;
00049     const realType *Acol = A, *Brow = B;
00050     for( i = 0; i < na * nc; ++i )
00051         C[i] = 0;
00052     for( k = 0; k < nb; ++k, Acol += na, Brow += nc )
00053     {
00054         realType* Ccol = C;
00055         for( j = 0; j < nc; ++j, Ccol += na )
00056             for( i = 0; i < na; ++i )
00057                 Ccol[i] += Acol[i] * Brow[j];
00058     }
00059 }
00060
00061 /*
00062 static void mxm_rr(const realType *A, unsigned na,
00063                    const realType *B, unsigned nb,
00064                          realType *C, unsigned nc)
00065 {
00066   unsigned i,j,k;
00067   realType *Ccol = C;
00068   const realType *Bcol = B;
00069   for(j=0;j<nc;++j,Ccol+=na,++Bcol) {
00070     const realType *Arow = A;
00071     for(i=0;i<na;++i,Arow+=nb) {
00072       const realType *Bkj = Bcol;
00073       Ccol[i]=0.0;
00074       for(k=0;k<nb;++k,Bkj+=nc)
00075         Ccol[i] += Arow[k] * *Bkj;
00076     }
00077   }
00078 }
00079 */
00080
00081 /*--------------------------------------------------------------------------
00082    Matrix-Vector Multiplication
00083
00084    mxv_f (y,ny,A,x,nx) :
00085       gives y = A x where A is ny x nx
00086       f := r | c   to indicate A is in row- or column- major format
00087   --------------------------------------------------------------------------*/
00088
00089 static void mxv_c( realType* y, unsigned ny, const realType* A, const realType* x, unsigned nx )
00090 {
00091     realType *yp = y, *y_end = y + ny;
00092     const realType* x_end = x + nx;
00093     realType xk           = *x;
00094     do
00095     {
00096         *yp++ = *A++ * xk;
00097     } while( yp != y_end );
00098     for( ++x; x != x_end; ++x )
00099     {
00100         xk = *x;
00101         yp = y;
00102         do
00103         {
00104             *yp++ += *A++ * xk;
00105         } while( yp != y_end );
00106     }
00107 }
00108
00109 static void mxv_r( realType* y, unsigned ny, const realType* A, const realType* x, unsigned nx )
00110 {
00111     realType* y_end       = y + ny;
00112     const realType* x_end = x + nx;
00113     do
00114     {
00115         const realType* xp = x;
00116         realType sum       = *A++ * *xp++;
00117         while( xp != x_end )
00118         {
00119             sum += *A++ * *xp++;
00120         }
00121         *y++ = sum;
00122     } while( y != y_end );
00123 }
00124
00125 /*--------------------------------------------------------------------------
00126    Vector-Vector Multiplication
00127
00128    inner (u,v,n) : inner product
00129   --------------------------------------------------------------------------*/
00130
00131 /* precondition: n>=1 */
00132 static realType inner( const realType* u, const realType* v, unsigned n )
00133 {
00134     const realType* u_end = u + n;
00135     realType sum          = *u++ * *v++;
00136     while( u != u_end )
00137     {
00138         sum += *u++ * *v++;
00139     }
00140     return sum;
00141 }
00142
00143 /*--------------------------------------------------------------------------
00144    1-,2-,3-d Tensor Application
00145
00146    the 3d case:
00147    tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2)
00148      gives v = [ R (x) S (x) T ] u
00149      where R is mr x nr, S is ms x ns, T is mt x nt,
00150        each in row- or column-major format according to f := r | c
00151      u is nr x ns x nt in column-major format (inner index is r)
00152      v is mr x ms x mt in column-major format (inner index is r)
00153   --------------------------------------------------------------------------*/
00154
00155 void tensor_c1( const realType* R, unsigned mr, unsigned nr, const realType* u, realType* v )
00156 {
00157     mxv_c( v, mr, R, u, nr );
00158 }
00159
00160 void tensor_r1( const realType* R, unsigned mr, unsigned nr, const realType* u, realType* v )
00161 {
00162     mxv_r( v, mr, R, u, nr );
00163 }
00164
00165 /* W holds mr*ns reals */
00166 void tensor_c2( const realType* R,
00167                 unsigned mr,
00168                 unsigned nr,
00169                 const realType* S,
00170                 unsigned ms,
00171                 unsigned ns,
00172                 const realType* u,
00173                 realType* v,
00174                 realType* W )
00175 {
00176     mxm_cc( R, mr, u, nr, W, ns );
00177     mxm_cr( W, mr, S, ns, v, ms );
00178 }
00179
00180 /* W holds mr*ns reals */
00181 void tensor_r2( const realType* R,
00182                 unsigned mr,
00183                 unsigned nr,
00184                 const realType* S,
00185                 unsigned ms,
00186                 unsigned ns,
00187                 const realType* u,
00188                 realType* v,
00189                 realType* W )
00190 {
00191     mxm_rc( R, mr, u, nr, W, ns );
00192     mxm_cc( W, mr, S, ns, v, ms );
00193 }
00194
00195 /* W holds mr*ns*nt reals,
00196    Z holds mr*ms*nt reals */
00197 void tensor_c3( const realType* R,
00198                 unsigned mr,
00199                 unsigned nr,
00200                 const realType* S,
00201                 unsigned ms,
00202                 unsigned ns,
00203                 const realType* T,
00204                 unsigned mt,
00205                 unsigned nt,
00206                 const realType* u,
00207                 realType* v,
00208                 realType* W,
00209                 realType* Z )
00210 {
00211     unsigned n, mrns = mr * ns, mrms = mr * ms;
00212     realType* Zp = Z;
00213     mxm_cc( R, mr, u, nr, W, ns * nt );
00214     for( n = 0; n < nt; ++n, W += mrns, Zp += mrms )
00215         mxm_cr( W, mr, S, ns, Zp, ms );
00216     mxm_cr( Z, mrms, T, nt, v, mt );
00217 }
00218
00219 /* W holds mr*ns*nt reals,
00220    Z holds mr*ms*nt reals */
00221 void tensor_r3( const realType* R,
00222                 unsigned mr,
00223                 unsigned nr,
00224                 const realType* S,
00225                 unsigned ms,
00226                 unsigned ns,
00227                 const realType* T,
00228                 unsigned mt,
00229                 unsigned nt,
00230                 const realType* u,
00231                 realType* v,
00232                 realType* W,
00233                 realType* Z )
00234 {
00235     unsigned n, mrns = mr * ns, mrms = mr * ms;
00236     realType* Zp = Z;
00237     mxm_rc( R, mr, u, nr, W, ns * nt );
00238     for( n = 0; n < nt; ++n, W += mrns, Zp += mrms )
00239         mxm_cc( W, mr, S, ns, Zp, ms );
00240     mxm_cc( Z, mrms, T, nt, v, mt );
00241 }
00242
00243 /*--------------------------------------------------------------------------
00244    1-,2-,3-d Tensor Application of Row Vectors (for Interpolation)
00245
00246    the 3d case:
00247    v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work)
00248    same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2):
00249      gives v = [ Jr (x) Js (x) Jt ] u
00250      where Jr, Js, Jt are row vectors (interpolation weights)
00251      u is nr x ns x nt in column-major format (inner index is r)
00252      v is a scalar
00253   --------------------------------------------------------------------------*/
00254
00255 realType tensor_i1( const realType* Jr, unsigned nr, const realType* u )
00256 {
00257     return inner( Jr, u, nr );
00258 }
00259
00260 /* work holds ns reals */
00261 realType tensor_i2( const realType* Jr,
00262                     unsigned nr,
00263                     const realType* Js,
00264                     unsigned ns,
00265                     const realType* u,
00266                     realType* work )
00267 {
00268     mxv_r( work, ns, u, Jr, nr );
00269     return inner( Js, work, ns );
00270 }
00271
00272 /* work holds ns*nt + nt reals */
00273 realType tensor_i3( const realType* Jr,
00274                     unsigned nr,
00275                     const realType* Js,
00276                     unsigned ns,
00277                     const realType* Jt,
00278                     unsigned nt,
00279                     const realType* u,
00280                     realType* work )
00281 {
00282     realType* work2 = work + nt;
00283     mxv_r( work2, ns * nt, u, Jr, nr );
00284     mxv_r( work, nt, work2, Js, ns );
00285     return inner( Jt, work, nt );
00286 }
00287
00288 /*--------------------------------------------------------------------------
00289    1-,2-,3-d Tensor Application of Row Vectors
00290              for simultaneous Interpolation and Gradient computation
00291
00292    the 3d case:
00293    v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work)
00294      gives v   = [ Jr (x) Js (x) Jt ] u
00295            g_0 = [ Dr (x) Js (x) Jt ] u
00296            g_1 = [ Jr (x) Ds (x) Jt ] u
00297            g_2 = [ Jr (x) Js (x) Dt ] u
00298      where Jr,Dr,Js,Ds,Jt,Dt are row vectors
00299        (interpolation & derivative weights)
00300      u is nr x ns x nt in column-major format (inner index is r)
00301      v is a scalar, g is an array of 3 reals
00302   --------------------------------------------------------------------------*/
00303
00304 realType tensor_ig1( const realType* Jr, const realType* Dr, unsigned nr, const realType* u, realType* g )
00305 {
00306     *g = inner( Dr, u, nr );
00307     return inner( Jr, u, nr );
00308 }
00309
00310 /* work holds 2*ns reals */
00311 realType tensor_ig2( const realType* Jr,
00312                      const realType* Dr,
00313                      unsigned nr,
00314                      const realType* Js,
00315                      const realType* Ds,
00316                      unsigned ns,
00317                      const realType* u,
00318                      realType* g,
00319                      realType* work )
00320 {
00321     realType *a = work, *ar = a + ns;
00322     mxv_r( a, ns, u, Jr, nr );
00323     mxv_r( ar, ns, u, Dr, nr );
00324     g[0] = inner( Js, ar, ns );
00325     g[1] = inner( Ds, a, ns );
00326     return inner( Js, a, ns );
00327 }
00328
00329 /* work holds 2*ns*nt + 3*ns reals */
00330 realType tensor_ig3( const realType* Jr,
00331                      const realType* Dr,
00332                      unsigned nr,
00333                      const realType* Js,
00334                      const realType* Ds,
00335                      unsigned ns,
00336                      const realType* Jt,
00337                      const realType* Dt,
00338                      unsigned nt,
00339                      const realType* u,
00340                      realType* g,
00341                      realType* work )
00342 {
00343     unsigned nsnt = ns * nt;
00344     realType *a = work, *ar = a + nsnt, *b = ar + nsnt, *br = b + ns, *bs = br + ns;
00345     mxv_r( a, nsnt, u, Jr, nr );
00346     mxv_r( ar, nsnt, u, Dr, nr );
00347     mxv_r( b, nt, a, Js, ns );
00348     mxv_r( br, nt, ar, Js, ns );
00349     mxv_r( bs, nt, a, Ds, ns );
00350     g[0] = inner( Jt, br, nt );
00351     g[1] = inner( Jt, bs, nt );
00352     g[2] = inner( Dt, b, nt );
00353     return inner( Jt, b, nt );
00354 }