MOAB: Mesh Oriented datABase  (version 5.3.1)
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, unsigned mr, unsigned nr, const realType* S, unsigned ms, unsigned ns,
00167                 const realType* u, realType* v, realType* W )
00168 {
00169     mxm_cc( R, mr, u, nr, W, ns );
00170     mxm_cr( W, mr, S, ns, v, ms );
00171 }
00172 
00173 /* W holds mr*ns reals */
00174 void tensor_r2( const realType* R, unsigned mr, unsigned nr, const realType* S, unsigned ms, unsigned ns,
00175                 const realType* u, realType* v, realType* W )
00176 {
00177     mxm_rc( R, mr, u, nr, W, ns );
00178     mxm_cc( W, mr, S, ns, v, ms );
00179 }
00180 
00181 /* W holds mr*ns*nt reals,
00182    Z holds mr*ms*nt reals */
00183 void tensor_c3( const realType* R, unsigned mr, unsigned nr, const realType* S, unsigned ms, unsigned ns,
00184                 const realType* T, unsigned mt, unsigned nt, const realType* u, realType* v, realType* W, realType* Z )
00185 {
00186     unsigned  n, mrns = mr * ns, mrms = mr * ms;
00187     realType* Zp = Z;
00188     mxm_cc( R, mr, u, nr, W, ns * nt );
00189     for( n = 0; n < nt; ++n, W += mrns, Zp += mrms )
00190         mxm_cr( W, mr, S, ns, Zp, ms );
00191     mxm_cr( Z, mrms, T, nt, v, mt );
00192 }
00193 
00194 /* W holds mr*ns*nt reals,
00195    Z holds mr*ms*nt reals */
00196 void tensor_r3( const realType* R, unsigned mr, unsigned nr, const realType* S, unsigned ms, unsigned ns,
00197                 const realType* T, unsigned mt, unsigned nt, const realType* u, realType* v, realType* W, realType* Z )
00198 {
00199     unsigned  n, mrns = mr * ns, mrms = mr * ms;
00200     realType* Zp = Z;
00201     mxm_rc( R, mr, u, nr, W, ns * nt );
00202     for( n = 0; n < nt; ++n, W += mrns, Zp += mrms )
00203         mxm_cc( W, mr, S, ns, Zp, ms );
00204     mxm_cc( Z, mrms, T, nt, v, mt );
00205 }
00206 
00207 /*--------------------------------------------------------------------------
00208    1-,2-,3-d Tensor Application of Row Vectors (for Interpolation)
00209 
00210    the 3d case:
00211    v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work)
00212    same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2):
00213      gives v = [ Jr (x) Js (x) Jt ] u
00214      where Jr, Js, Jt are row vectors (interpolation weights)
00215      u is nr x ns x nt in column-major format (inner index is r)
00216      v is a scalar
00217   --------------------------------------------------------------------------*/
00218 
00219 realType tensor_i1( const realType* Jr, unsigned nr, const realType* u )
00220 {
00221     return inner( Jr, u, nr );
00222 }
00223 
00224 /* work holds ns reals */
00225 realType tensor_i2( const realType* Jr, unsigned nr, const realType* Js, unsigned ns, const realType* u,
00226                     realType* work )
00227 {
00228     mxv_r( work, ns, u, Jr, nr );
00229     return inner( Js, work, ns );
00230 }
00231 
00232 /* work holds ns*nt + nt reals */
00233 realType tensor_i3( const realType* Jr, unsigned nr, const realType* Js, unsigned ns, const realType* Jt, unsigned nt,
00234                     const realType* u, realType* work )
00235 {
00236     realType* work2 = work + nt;
00237     mxv_r( work2, ns * nt, u, Jr, nr );
00238     mxv_r( work, nt, work2, Js, ns );
00239     return inner( Jt, work, nt );
00240 }
00241 
00242 /*--------------------------------------------------------------------------
00243    1-,2-,3-d Tensor Application of Row Vectors
00244              for simultaneous Interpolation and Gradient computation
00245 
00246    the 3d case:
00247    v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work)
00248      gives v   = [ Jr (x) Js (x) Jt ] u
00249            g_0 = [ Dr (x) Js (x) Jt ] u
00250            g_1 = [ Jr (x) Ds (x) Jt ] u
00251            g_2 = [ Jr (x) Js (x) Dt ] u
00252      where Jr,Dr,Js,Ds,Jt,Dt are row vectors
00253        (interpolation & derivative weights)
00254      u is nr x ns x nt in column-major format (inner index is r)
00255      v is a scalar, g is an array of 3 reals
00256   --------------------------------------------------------------------------*/
00257 
00258 realType tensor_ig1( const realType* Jr, const realType* Dr, unsigned nr, const realType* u, realType* g )
00259 {
00260     *g = inner( Dr, u, nr );
00261     return inner( Jr, u, nr );
00262 }
00263 
00264 /* work holds 2*ns reals */
00265 realType tensor_ig2( const realType* Jr, const realType* Dr, unsigned nr, const realType* Js, const realType* Ds,
00266                      unsigned ns, const realType* u, realType* g, realType* work )
00267 {
00268     realType *a = work, *ar = a + ns;
00269     mxv_r( a, ns, u, Jr, nr );
00270     mxv_r( ar, ns, u, Dr, nr );
00271     g[ 0 ] = inner( Js, ar, ns );
00272     g[ 1 ] = inner( Ds, a, ns );
00273     return inner( Js, a, ns );
00274 }
00275 
00276 /* work holds 2*ns*nt + 3*ns reals */
00277 realType tensor_ig3( const realType* Jr, const realType* Dr, unsigned nr, const realType* Js, const realType* Ds,
00278                      unsigned ns, const realType* Jt, const realType* Dt, unsigned nt, const realType* u, realType* g,
00279                      realType* work )
00280 {
00281     unsigned  nsnt = ns * nt;
00282     realType *a = work, *ar = a + nsnt, *b = ar + nsnt, *br = b + ns, *bs = br + ns;
00283     mxv_r( a, nsnt, u, Jr, nr );
00284     mxv_r( ar, nsnt, u, Dr, nr );
00285     mxv_r( b, nt, a, Js, ns );
00286     mxv_r( br, nt, ar, Js, ns );
00287     mxv_r( bs, nt, a, Ds, ns );
00288     g[ 0 ] = inner( Jt, br, nt );
00289     g[ 1 ] = inner( Jt, bs, nt );
00290     g[ 2 ] = inner( Dt, b, nt );
00291     return inner( Jt, b, nt );
00292 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines