MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 }