/[escript]/trunk/paso/src/SparseMatrix_MatrixMatrixTranspose.cpp
ViewVC logotype

Diff of /trunk/paso/src/SparseMatrix_MatrixMatrixTranspose.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4828 by caltinay, Tue Apr 1 03:50:23 2014 UTC revision 4829 by caltinay, Thu Apr 3 04:02:53 2014 UTC
# Line 15  Line 15 
15  *****************************************************************************/  *****************************************************************************/
16    
17    
18  /************************************************************************************  /****************************************************************************
19    
20   Paso: Sparse matrix product (for efficiency, use the transpose   Paso: Sparse matrix product (for efficiency, use the transpose
21         of Matrix B when B^T is available)         of Matrix B when B^T is available)
22    
23  ************************************************************************************  *****************************************************************************
24    
25     Author: l.gao@uq.edu.au     Author: l.gao@uq.edu.au
26    
27  ************************************************************************************/  *****************************************************************************/
28    
29  #include "SparseMatrix.h"  #include "SparseMatrix.h"
30  #include "Paso.h"  #include "Paso.h"
# Line 32  Line 32 
32  namespace paso {  namespace paso {
33    
34  // forward declarations  // forward declarations
35  void SparseMatrix_MatrixMatrixTranspose_DD(SparseMatrix *C,  void SparseMatrix_MatrixMatrixTranspose_DD(SparseMatrix_ptr C,
36          const SparseMatrix* A, const SparseMatrix* B, const SparseMatrix* T);                                             const_SparseMatrix_ptr A,
37  void SparseMatrix_MatrixMatrixTranspose_DB(SparseMatrix *C,                                             const_SparseMatrix_ptr B,
38          const SparseMatrix* A, const SparseMatrix* B, const SparseMatrix* T);                                             const_SparseMatrix_ptr T);
39  void SparseMatrix_MatrixMatrixTranspose_BD(SparseMatrix *C,  void SparseMatrix_MatrixMatrixTranspose_DB(SparseMatrix_ptr C,
40          const SparseMatrix* A, const SparseMatrix* B, const SparseMatrix* T);                                             const_SparseMatrix_ptr A,
41  void SparseMatrix_MatrixMatrixTranspose_BB(SparseMatrix *C,                                             const_SparseMatrix_ptr B,
42          const SparseMatrix* A, const SparseMatrix* B, const SparseMatrix* T);                                             const_SparseMatrix_ptr T);
43    void SparseMatrix_MatrixMatrixTranspose_BD(SparseMatrix_ptr C,
44  SparseMatrix* SparseMatrix_MatrixMatrixTranspose(const SparseMatrix* A, const SparseMatrix* B, const SparseMatrix* T)                                             const_SparseMatrix_ptr A,
45                                               const_SparseMatrix_ptr B,
46                                               const_SparseMatrix_ptr T);
47    void SparseMatrix_MatrixMatrixTranspose_BB(SparseMatrix_ptr C,
48                                               const_SparseMatrix_ptr A,
49                                               const_SparseMatrix_ptr B,
50                                               const_SparseMatrix_ptr T);
51    
52    SparseMatrix_ptr SparseMatrix_MatrixMatrixTranspose(const_SparseMatrix_ptr A,
53                                                        const_SparseMatrix_ptr B,
54                                                        const_SparseMatrix_ptr T)
55  {  {
56     SparseMatrixType C_type;      SparseMatrixType C_type;
57          SparseMatrix_ptr out;
    SparseMatrix *out=NULL;  
58    
59     if ( !  ( (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) || (A->type & MATRIX_FORMAT_DEFAULT) || (MATRIX_FORMAT_BLK1 & A->type ) )  ) {      if ( !  ( (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) || (A->type & MATRIX_FORMAT_DEFAULT) || (MATRIX_FORMAT_BLK1 & A->type ) )  ) {
60        Esys_setError(TYPE_ERROR,"SparseMatrix_MatrixMatrix: Unsupported matrix format of A.");          Esys_setError(TYPE_ERROR,"SparseMatrix_MatrixMatrix: Unsupported matrix format of A.");
61        return NULL;          return out;
62     }      }
63     if ( !  ( (B->type & MATRIX_FORMAT_DIAGONAL_BLOCK) || (B->type & MATRIX_FORMAT_DEFAULT) || (MATRIX_FORMAT_BLK1 & B->type ) ) ) {      if ( !  ( (B->type & MATRIX_FORMAT_DIAGONAL_BLOCK) || (B->type & MATRIX_FORMAT_DEFAULT) || (MATRIX_FORMAT_BLK1 & B->type ) ) ) {
64        Esys_setError(TYPE_ERROR,"SparseMatrix_MatrixMatrix: Unsupported matrix format of B.");          Esys_setError(TYPE_ERROR,"SparseMatrix_MatrixMatrix: Unsupported matrix format of B.");
65        return NULL;          return out;
66     }      }
67     if (! (A->col_block_size == B->row_block_size) ) {      if (! (A->col_block_size == B->row_block_size) ) {
68        Esys_setError(TYPE_ERROR,"SparseMatrix_MatrixMatrix: Column block size of A and row block size of B must match.");          Esys_setError(TYPE_ERROR,"SparseMatrix_MatrixMatrix: Column block size of A and row block size of B must match.");
69        return NULL;          return out;
70     }      }
71     if (! (A->numCols == B->numRows) ) {      if (! (A->numCols == B->numRows) ) {
72        Esys_setError(TYPE_ERROR,"SparseMatrix_MatrixMatrix: number of columns of A and number of rows of B must match.");          Esys_setError(TYPE_ERROR,"SparseMatrix_MatrixMatrix: number of columns of A and number of rows of B must match.");
73        return NULL;          return out;
74     }      }
75      
76          if ( (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) && (B->type & MATRIX_FORMAT_DIAGONAL_BLOCK) ) {
77     if ( (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) && (B->type & MATRIX_FORMAT_DIAGONAL_BLOCK) ) {          C_type=MATRIX_FORMAT_DIAGONAL_BLOCK;
78        C_type=MATRIX_FORMAT_DIAGONAL_BLOCK;      } else {
79     } else {          C_type=MATRIX_FORMAT_DEFAULT;
80        C_type=MATRIX_FORMAT_DEFAULT;      }
    }  
81    
82     Pattern_ptr outpattern(A->pattern->multiply(MATRIX_FORMAT_DEFAULT,B->pattern));      Pattern_ptr outpattern(A->pattern->multiply(MATRIX_FORMAT_DEFAULT, B->pattern));
83        
84     if (Esys_noError()) {      if (Esys_noError()) {
85        out=SparseMatrix_alloc(C_type, outpattern, A->row_block_size, B->col_block_size, FALSE);          out.reset(new SparseMatrix(C_type, outpattern, A->row_block_size, B->col_block_size, false));
86     }      }
87    
88     if (Esys_noError()) {      if (Esys_noError()) {
89        if ( (A->row_block_size == 1) && (B->col_block_size ==1 ) && (A->col_block_size ==1) ) {          if ( (A->row_block_size == 1) && (B->col_block_size ==1 ) && (A->col_block_size ==1) ) {
90       SparseMatrix_MatrixMatrixTranspose_DD(out, A,B,T);              SparseMatrix_MatrixMatrixTranspose_DD(out, A, B, T);
91        } else {          } else {
92       if (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) {              if (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) {
93          if (B->type & MATRIX_FORMAT_DIAGONAL_BLOCK) {                  if (B->type & MATRIX_FORMAT_DIAGONAL_BLOCK) {
94             SparseMatrix_MatrixMatrixTranspose_DD(out, A,B,T);                      SparseMatrix_MatrixMatrixTranspose_DD(out, A, B, T);
95          } else {                  } else {
96             SparseMatrix_MatrixMatrixTranspose_DB(out, A,B,T);                      SparseMatrix_MatrixMatrixTranspose_DB(out, A, B, T);
97          }                  }
98       } else {              } else {
99          if (B->type & MATRIX_FORMAT_DIAGONAL_BLOCK) {                  if (B->type & MATRIX_FORMAT_DIAGONAL_BLOCK) {
100             SparseMatrix_MatrixMatrixTranspose_BD(out, A,B,T);                      SparseMatrix_MatrixMatrixTranspose_BD(out, A, B, T);
101          } else {                  } else {
102             SparseMatrix_MatrixMatrixTranspose_BB(out, A,B,T);                      SparseMatrix_MatrixMatrixTranspose_BB(out, A, B, T);
103          }                  }
104       }              }
105        }          }
106        return out;      } else {
107     } else {          out.reset();
108        SparseMatrix_free(out);      }
109        return NULL;      return out;
    }  
110  }  }
111    
112  /* not good for block size 1 */  /* not good for block size 1 */
113  void SparseMatrix_MatrixMatrixTranspose_BB(SparseMatrix* C, const SparseMatrix* A, const SparseMatrix* B, const SparseMatrix* T)  void SparseMatrix_MatrixMatrixTranspose_BB(SparseMatrix_ptr C, const_SparseMatrix_ptr A, const_SparseMatrix_ptr B, const_SparseMatrix_ptr T)
114  {  {
115     const dim_t n = C->numRows;     const dim_t n = C->numRows;
116     const dim_t row_block_size = C->row_block_size;     const dim_t row_block_size = C->row_block_size;
# Line 120  void SparseMatrix_MatrixMatrixTranspose_ Line 127  void SparseMatrix_MatrixMatrixTranspose_
127     if ( (row_block_size == 2) && (col_block_size ==2 ) && (A_col_block_size == 2) ) {     if ( (row_block_size == 2) && (col_block_size ==2 ) && (A_col_block_size == 2) ) {
128        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_01, C_ij_11)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_01, C_ij_11)
129        {        {
130       #pragma omp for schedule(static)           #pragma omp for schedule(static)
131       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
132          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
133             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
134                            
135             C_ij_00=0;                 C_ij_00=0;
136             C_ij_10=0;                 C_ij_10=0;
137             C_ij_01=0;                 C_ij_01=0;
138             C_ij_11=0;                 C_ij_11=0;
139    
140             C_ij=&(C->val[ij_ptrC*4]);                 C_ij=&(C->val[ij_ptrC*4]);
141    
142                 ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
143                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
# Line 161  void SparseMatrix_MatrixMatrixTranspose_ Line 168  void SparseMatrix_MatrixMatrixTranspose_
168                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
169                   }                   }
170                 }                 }
171             C_ij[0+2*0]=C_ij_00;                 C_ij[0+2*0]=C_ij_00;
172             C_ij[1+2*0]=C_ij_10;                 C_ij[1+2*0]=C_ij_10;
173             C_ij[0+2*1]=C_ij_01;                 C_ij[0+2*1]=C_ij_01;
174             C_ij[1+2*1]=C_ij_11;                 C_ij[1+2*1]=C_ij_11;
175    
176          }              }
177       }           }
178        } /* end of parallel region */        } /* end of parallel region */
179                
180     } else if ( (row_block_size == 3) && (col_block_size ==3 ) && (A_col_block_size == 3)  ){     } else if ( (row_block_size == 3) && (col_block_size ==3 ) && (A_col_block_size == 3)  ){
181        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_01, C_ij_11, C_ij_21, C_ij_02, C_ij_12, C_ij_22)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_01, C_ij_11, C_ij_21, C_ij_02, C_ij_12, C_ij_22)
182        {        {
183       #pragma omp for schedule(static)           #pragma omp for schedule(static)
184       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
185          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
186             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
187                            
188             C_ij_00=0;                 C_ij_00=0;
189             C_ij_10=0;                 C_ij_10=0;
190             C_ij_20=0;                 C_ij_20=0;
191             C_ij_01=0;                 C_ij_01=0;
192             C_ij_11=0;                 C_ij_11=0;
193             C_ij_21=0;                 C_ij_21=0;
194             C_ij_02=0;                 C_ij_02=0;
195             C_ij_12=0;                 C_ij_12=0;
196             C_ij_22=0;                 C_ij_22=0;
197                            
198             C_ij=&(C->val[ij_ptrC*9]);                 C_ij=&(C->val[ij_ptrC*9]);
199    
200             ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
201                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
202                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
203                 kjb = T->pattern->ptr[j+1];                 kjb = T->pattern->ptr[j+1];
# Line 206  void SparseMatrix_MatrixMatrixTranspose_ Line 213  void SparseMatrix_MatrixMatrixTranspose_
213                     if (kj_ptrB >= kjb) break;                     if (kj_ptrB >= kjb) break;
214                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
215                   } else {                   } else {
216              A_ik=&(A->val[ik_ptrA*9]);                          A_ik=&(A->val[ik_ptrA*9]);
217              B_kj=&(T->val[kj_ptrB*9]);                          B_kj=&(T->val[kj_ptrB*9]);
218                                        
219              C_ij_00 +=A_ik[0+3*0]*B_kj[0+3*0]                          C_ij_00 +=A_ik[0+3*0]*B_kj[0+3*0]
220                       +A_ik[0+3*1]*B_kj[1+3*0]                                   +A_ik[0+3*1]*B_kj[1+3*0]
221                       +A_ik[0+3*2]*B_kj[2+3*0];                                   +A_ik[0+3*2]*B_kj[2+3*0];
222              C_ij_10 +=A_ik[1+3*0]*B_kj[0+3*0]                          C_ij_10 +=A_ik[1+3*0]*B_kj[0+3*0]
223                       +A_ik[1+3*1]*B_kj[1+3*0]                                   +A_ik[1+3*1]*B_kj[1+3*0]
224                       +A_ik[1+3*2]*B_kj[2+3*0];                                   +A_ik[1+3*2]*B_kj[2+3*0];
225              C_ij_20 +=A_ik[2+3*0]*B_kj[0+3*0]                          C_ij_20 +=A_ik[2+3*0]*B_kj[0+3*0]
226                   +A_ik[2+3*1]*B_kj[1+3*0]                                   +A_ik[2+3*1]*B_kj[1+3*0]
227                   +A_ik[2+3*2]*B_kj[2+3*0];                                   +A_ik[2+3*2]*B_kj[2+3*0];
228                                                    
229              C_ij_01 +=A_ik[0+3*0]*B_kj[0+3*1]                          C_ij_01 +=A_ik[0+3*0]*B_kj[0+3*1]
230                       +A_ik[0+3*1]*B_kj[1+3*1]                                   +A_ik[0+3*1]*B_kj[1+3*1]
231                       +A_ik[0+3*2]*B_kj[2+3*1];                                   +A_ik[0+3*2]*B_kj[2+3*1];
232              C_ij_11 +=A_ik[1+3*0]*B_kj[0+3*1]                          C_ij_11 +=A_ik[1+3*0]*B_kj[0+3*1]
233                       +A_ik[1+3*1]*B_kj[1+3*1]                                   +A_ik[1+3*1]*B_kj[1+3*1]
234                       +A_ik[1+3*2]*B_kj[2+3*1];                                   +A_ik[1+3*2]*B_kj[2+3*1];
235              C_ij_21 +=A_ik[2+3*0]*B_kj[0+3*1]                          C_ij_21 +=A_ik[2+3*0]*B_kj[0+3*1]
236                   +A_ik[2+3*1]*B_kj[1+3*1]                                   +A_ik[2+3*1]*B_kj[1+3*1]
237                   +A_ik[2+3*2]*B_kj[2+3*1];                                   +A_ik[2+3*2]*B_kj[2+3*1];
238    
239              C_ij_01 +=A_ik[0+3*0]*B_kj[0+3*2]                          C_ij_01 +=A_ik[0+3*0]*B_kj[0+3*2]
240                       +A_ik[0+3*1]*B_kj[1+3*2]                                   +A_ik[0+3*1]*B_kj[1+3*2]
241                       +A_ik[0+3*2]*B_kj[2+3*2];                                   +A_ik[0+3*2]*B_kj[2+3*2];
242              C_ij_11 +=A_ik[1+3*0]*B_kj[0+3*2]                          C_ij_11 +=A_ik[1+3*0]*B_kj[0+3*2]
243                       +A_ik[1+3*1]*B_kj[1+3*2]                                   +A_ik[1+3*1]*B_kj[1+3*2]
244                       +A_ik[1+3*2]*B_kj[2+3*2];                                   +A_ik[1+3*2]*B_kj[2+3*2];
245              C_ij_21 +=A_ik[2+3*0]*B_kj[0+3*2]                          C_ij_21 +=A_ik[2+3*0]*B_kj[0+3*2]
246                   +A_ik[2+3*1]*B_kj[1+3*2]                                   +A_ik[2+3*1]*B_kj[1+3*2]
247                   +A_ik[2+3*2]*B_kj[2+3*2];                                   +A_ik[2+3*2]*B_kj[2+3*2];
248              ik_ptrA ++;                          ik_ptrA ++;
249              kj_ptrB ++;                          kj_ptrB ++;
250              if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                          if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
251              kA=A->pattern->index[ik_ptrA];                          kA=A->pattern->index[ik_ptrA];
252              kB=T->pattern->index[kj_ptrB];                          kB=T->pattern->index[kj_ptrB];
253           }                   }
254             }                 }
255             C_ij[0+3*0]=C_ij_00;                 C_ij[0+3*0]=C_ij_00;
256             C_ij[1+3*0]=C_ij_10;                 C_ij[1+3*0]=C_ij_10;
257             C_ij[2+3*0]=C_ij_20;                 C_ij[2+3*0]=C_ij_20;
258             C_ij[0+3*1]=C_ij_01;                 C_ij[0+3*1]=C_ij_01;
259             C_ij[1+3*1]=C_ij_11;                 C_ij[1+3*1]=C_ij_11;
260             C_ij[2+3*1]=C_ij_21;                 C_ij[2+3*1]=C_ij_21;
261             C_ij[0+3*2]=C_ij_02;                 C_ij[0+3*2]=C_ij_02;
262             C_ij[1+3*2]=C_ij_12;                 C_ij[1+3*2]=C_ij_12;
263             C_ij[2+3*2]=C_ij_22;                 C_ij[2+3*2]=C_ij_22;
264          }              }
265       }           }
266        } /* end of parallel region */        } /* end of parallel region */
267     } else if ( (row_block_size == 4) && (col_block_size ==4 ) && (A_col_block_size == 4)  ){     } else if ( (row_block_size == 4) && (col_block_size ==4 ) && (A_col_block_size == 4)  ){
268        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_30, C_ij_01, C_ij_11, C_ij_21, C_ij_31, C_ij_02, C_ij_12, C_ij_22, C_ij_32, C_ij_03, C_ij_13, C_ij_23, C_ij_33)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_30, C_ij_01, C_ij_11, C_ij_21, C_ij_31, C_ij_02, C_ij_12, C_ij_22, C_ij_32, C_ij_03, C_ij_13, C_ij_23, C_ij_33)
269        {        {
270       #pragma omp for schedule(static)           #pragma omp for schedule(static)
271       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
272          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
273             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
274                            
275             C_ij_00=0;                 C_ij_00=0;
276             C_ij_10=0;                 C_ij_10=0;
277             C_ij_20=0;                 C_ij_20=0;
278             C_ij_30=0;                 C_ij_30=0;
279             C_ij_01=0;                 C_ij_01=0;
280             C_ij_11=0;                 C_ij_11=0;
281             C_ij_21=0;                 C_ij_21=0;
282             C_ij_31=0;                 C_ij_31=0;
283             C_ij_02=0;                 C_ij_02=0;
284             C_ij_12=0;                 C_ij_12=0;
285             C_ij_22=0;                 C_ij_22=0;
286             C_ij_32=0;                 C_ij_32=0;
287             C_ij_03=0;                 C_ij_03=0;
288             C_ij_13=0;                 C_ij_13=0;
289             C_ij_23=0;                 C_ij_23=0;
290             C_ij_33=0;                           C_ij_33=0;              
291                            
292             C_ij=&(C->val[ij_ptrC*16]);                 C_ij=&(C->val[ij_ptrC*16]);
293    
294                 ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
295                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
# Line 300  void SparseMatrix_MatrixMatrixTranspose_ Line 307  void SparseMatrix_MatrixMatrixTranspose_
307                     if (kj_ptrB >= kjb) break;                     if (kj_ptrB >= kjb) break;
308                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
309                   } else {                   } else {
310                 A_ik=&(A->val[ik_ptrA*16]);                             A_ik=&(A->val[ik_ptrA*16]);
311                 B_kj=&(T->val[kj_ptrB*16]);                             B_kj=&(T->val[kj_ptrB*16]);
312                                            
313                 C_ij_00 +=A_ik[0+4*0]*B_kj[0+4*0]                             C_ij_00 +=A_ik[0+4*0]*B_kj[0+4*0]
314                      +A_ik[0+4*1]*B_kj[1+4*0]                                      +A_ik[0+4*1]*B_kj[1+4*0]
315                      +A_ik[0+4*2]*B_kj[2+4*0]                                      +A_ik[0+4*2]*B_kj[2+4*0]
316                      +A_ik[0+4*3]*B_kj[3+4*0];                                      +A_ik[0+4*3]*B_kj[3+4*0];
317                 C_ij_10 +=A_ik[1+4*0]*B_kj[0+4*0]                             C_ij_10 +=A_ik[1+4*0]*B_kj[0+4*0]
318                      +A_ik[1+4*1]*B_kj[1+4*0]                                      +A_ik[1+4*1]*B_kj[1+4*0]
319                      +A_ik[1+4*2]*B_kj[2+4*0]                                      +A_ik[1+4*2]*B_kj[2+4*0]
320                      +A_ik[1+4*3]*B_kj[3+4*0];                                      +A_ik[1+4*3]*B_kj[3+4*0];
321                 C_ij_20 +=A_ik[2+4*0]*B_kj[0+4*0]                             C_ij_20 +=A_ik[2+4*0]*B_kj[0+4*0]
322                      +A_ik[2+4*1]*B_kj[1+4*0]                                      +A_ik[2+4*1]*B_kj[1+4*0]
323                      +A_ik[2+4*2]*B_kj[2+4*0]                                      +A_ik[2+4*2]*B_kj[2+4*0]
324                      +A_ik[2+4*3]*B_kj[3+4*0];                                      +A_ik[2+4*3]*B_kj[3+4*0];
325                 C_ij_30 +=A_ik[3+4*0]*B_kj[0+4*0]                             C_ij_30 +=A_ik[3+4*0]*B_kj[0+4*0]
326                      +A_ik[3+4*1]*B_kj[1+4*0]                                      +A_ik[3+4*1]*B_kj[1+4*0]
327                      +A_ik[3+4*2]*B_kj[2+4*0]                                      +A_ik[3+4*2]*B_kj[2+4*0]
328                      +A_ik[3+4*3]*B_kj[3+4*0];                                      +A_ik[3+4*3]*B_kj[3+4*0];
329                                                            
330                 C_ij_01 +=A_ik[0+4*0]*B_kj[0+4*1]                             C_ij_01 +=A_ik[0+4*0]*B_kj[0+4*1]
331                      +A_ik[0+4*1]*B_kj[1+4*1]                                      +A_ik[0+4*1]*B_kj[1+4*1]
332                      +A_ik[0+4*2]*B_kj[2+4*1]                                      +A_ik[0+4*2]*B_kj[2+4*1]
333                      +A_ik[0+4*3]*B_kj[3+4*1];                                      +A_ik[0+4*3]*B_kj[3+4*1];
334                 C_ij_11 +=A_ik[1+4*0]*B_kj[0+4*1]                             C_ij_11 +=A_ik[1+4*0]*B_kj[0+4*1]
335                      +A_ik[1+4*1]*B_kj[1+4*1]                                      +A_ik[1+4*1]*B_kj[1+4*1]
336                      +A_ik[1+4*2]*B_kj[2+4*1]                                      +A_ik[1+4*2]*B_kj[2+4*1]
337                      +A_ik[1+4*3]*B_kj[3+4*1];                                      +A_ik[1+4*3]*B_kj[3+4*1];
338                 C_ij_21 +=A_ik[2+4*0]*B_kj[0+4*1]                             C_ij_21 +=A_ik[2+4*0]*B_kj[0+4*1]
339                      +A_ik[2+4*1]*B_kj[1+4*1]                                      +A_ik[2+4*1]*B_kj[1+4*1]
340                      +A_ik[2+4*2]*B_kj[2+4*1]                                      +A_ik[2+4*2]*B_kj[2+4*1]
341                      +A_ik[2+4*3]*B_kj[3+4*1];                                      +A_ik[2+4*3]*B_kj[3+4*1];
342                 C_ij_31 +=A_ik[3+4*0]*B_kj[0+4*1]                             C_ij_31 +=A_ik[3+4*0]*B_kj[0+4*1]
343                      +A_ik[3+4*1]*B_kj[1+4*1]                                      +A_ik[3+4*1]*B_kj[1+4*1]
344                      +A_ik[3+4*2]*B_kj[2+4*1]                                      +A_ik[3+4*2]*B_kj[2+4*1]
345                      +A_ik[3+4*3]*B_kj[3+4*1];                                      +A_ik[3+4*3]*B_kj[3+4*1];
346    
347                 C_ij_02 +=A_ik[0+4*0]*B_kj[0+4*2]                             C_ij_02 +=A_ik[0+4*0]*B_kj[0+4*2]
348                      +A_ik[0+4*1]*B_kj[1+4*2]                                      +A_ik[0+4*1]*B_kj[1+4*2]
349                      +A_ik[0+4*2]*B_kj[2+4*2]                                      +A_ik[0+4*2]*B_kj[2+4*2]
350                      +A_ik[0+4*3]*B_kj[3+4*2];                                      +A_ik[0+4*3]*B_kj[3+4*2];
351                 C_ij_12 +=A_ik[1+4*0]*B_kj[0+4*2]                             C_ij_12 +=A_ik[1+4*0]*B_kj[0+4*2]
352                      +A_ik[1+4*1]*B_kj[1+4*2]                                      +A_ik[1+4*1]*B_kj[1+4*2]
353                      +A_ik[1+4*2]*B_kj[2+4*2]                                      +A_ik[1+4*2]*B_kj[2+4*2]
354                      +A_ik[1+4*3]*B_kj[3+4*2];                                      +A_ik[1+4*3]*B_kj[3+4*2];
355                 C_ij_22 +=A_ik[2+4*0]*B_kj[0+4*2]                             C_ij_22 +=A_ik[2+4*0]*B_kj[0+4*2]
356                      +A_ik[2+4*1]*B_kj[1+4*2]                                      +A_ik[2+4*1]*B_kj[1+4*2]
357                      +A_ik[2+4*2]*B_kj[2+4*2]                                      +A_ik[2+4*2]*B_kj[2+4*2]
358                      +A_ik[2+4*3]*B_kj[3+4*2];                                      +A_ik[2+4*3]*B_kj[3+4*2];
359                 C_ij_32 +=A_ik[3+4*0]*B_kj[0+4*2]                             C_ij_32 +=A_ik[3+4*0]*B_kj[0+4*2]
360                      +A_ik[3+4*1]*B_kj[1+4*2]                                      +A_ik[3+4*1]*B_kj[1+4*2]
361                      +A_ik[3+4*2]*B_kj[2+4*2]                                      +A_ik[3+4*2]*B_kj[2+4*2]
362                      +A_ik[3+4*3]*B_kj[3+4*2];                                      +A_ik[3+4*3]*B_kj[3+4*2];
363    
364                 C_ij_03 +=A_ik[0+4*0]*B_kj[0+4*3]                             C_ij_03 +=A_ik[0+4*0]*B_kj[0+4*3]
365                      +A_ik[0+4*1]*B_kj[1+4*3]                                      +A_ik[0+4*1]*B_kj[1+4*3]
366                      +A_ik[0+4*2]*B_kj[2+4*3]                                      +A_ik[0+4*2]*B_kj[2+4*3]
367                      +A_ik[0+4*3]*B_kj[3+4*3];                                      +A_ik[0+4*3]*B_kj[3+4*3];
368                 C_ij_13 +=A_ik[1+4*0]*B_kj[0+4*3]                             C_ij_13 +=A_ik[1+4*0]*B_kj[0+4*3]
369                      +A_ik[1+4*1]*B_kj[1+4*3]                                      +A_ik[1+4*1]*B_kj[1+4*3]
370                      +A_ik[1+4*2]*B_kj[2+4*3]                                      +A_ik[1+4*2]*B_kj[2+4*3]
371                      +A_ik[1+4*3]*B_kj[3+4*3];                                      +A_ik[1+4*3]*B_kj[3+4*3];
372                 C_ij_23 +=A_ik[2+4*0]*B_kj[0+4*3]                             C_ij_23 +=A_ik[2+4*0]*B_kj[0+4*3]
373                      +A_ik[2+4*1]*B_kj[1+4*3]                                      +A_ik[2+4*1]*B_kj[1+4*3]
374                      +A_ik[2+4*2]*B_kj[2+4*3]                                      +A_ik[2+4*2]*B_kj[2+4*3]
375                      +A_ik[2+4*3]*B_kj[3+4*3];                                      +A_ik[2+4*3]*B_kj[3+4*3];
376                 C_ij_33 +=A_ik[3+4*0]*B_kj[0+4*3]                             C_ij_33 +=A_ik[3+4*0]*B_kj[0+4*3]
377                      +A_ik[3+4*1]*B_kj[1+4*3]                                      +A_ik[3+4*1]*B_kj[1+4*3]
378                      +A_ik[3+4*2]*B_kj[2+4*3]                                      +A_ik[3+4*2]*B_kj[2+4*3]
379                      +A_ik[3+4*3]*B_kj[3+4*3];                                      +A_ik[3+4*3]*B_kj[3+4*3];
380                 ik_ptrA ++;                             ik_ptrA ++;
381                 kj_ptrB ++;                             kj_ptrB ++;
382                 if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                             if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
383                 kA=A->pattern->index[ik_ptrA];                             kA=A->pattern->index[ik_ptrA];
384                 kB=T->pattern->index[kj_ptrB];                             kB=T->pattern->index[kj_ptrB];
385            }                    }
386          }                  }
387          C_ij[0+4*0]=C_ij_00;                  C_ij[0+4*0]=C_ij_00;
388          C_ij[1+4*0]=C_ij_10;                  C_ij[1+4*0]=C_ij_10;
389          C_ij[2+4*0]=C_ij_20;                  C_ij[2+4*0]=C_ij_20;
390          C_ij[3+4*0]=C_ij_30;                  C_ij[3+4*0]=C_ij_30;
391          C_ij[0+4*1]=C_ij_01;                  C_ij[0+4*1]=C_ij_01;
392          C_ij[1+4*1]=C_ij_11;                  C_ij[1+4*1]=C_ij_11;
393          C_ij[2+4*1]=C_ij_21;                  C_ij[2+4*1]=C_ij_21;
394          C_ij[3+4*1]=C_ij_31;                  C_ij[3+4*1]=C_ij_31;
395          C_ij[0+4*2]=C_ij_02;                  C_ij[0+4*2]=C_ij_02;
396          C_ij[1+4*2]=C_ij_12;                  C_ij[1+4*2]=C_ij_12;
397          C_ij[2+4*2]=C_ij_22;                  C_ij[2+4*2]=C_ij_22;
398          C_ij[3+4*2]=C_ij_32;                  C_ij[3+4*2]=C_ij_32;
399          C_ij[0+4*3]=C_ij_03;                  C_ij[0+4*3]=C_ij_03;
400          C_ij[1+4*3]=C_ij_13;                  C_ij[1+4*3]=C_ij_13;
401          C_ij[2+4*3]=C_ij_23;                  C_ij[2+4*3]=C_ij_23;
402          C_ij[3+4*3]=C_ij_33;                            C_ij[3+4*3]=C_ij_33;              
403          }              }
404       }           }
405        } /* end of parallel region */        } /* end of parallel region */
406                    
407     } else {     } else {
408        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, rtmp, A_ik,B_kj )        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, rtmp, A_ik,B_kj )
409        {        {
410       #pragma omp for schedule(static)           #pragma omp for schedule(static)
411       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
412          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
413             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
414             C_ij=&(C->val[ij_ptrC*C_block_size]);                 C_ij=&(C->val[ij_ptrC*C_block_size]);
415             for (ib=0; ib<C_block_size; ++ib)  C_ij[ib]=0;                 for (ib=0; ib<C_block_size; ++ib)  C_ij[ib]=0;
416                 ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
417                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
418                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
# Line 424  void SparseMatrix_MatrixMatrixTranspose_ Line 431  void SparseMatrix_MatrixMatrixTranspose_
431                   } else {                   } else {
432                     A_ik=&(A->val[ik_ptrA*A_block_size]);                     A_ik=&(A->val[ik_ptrA*A_block_size]);
433                     B_kj=&(T->val[kj_ptrB*B_block_size]);                     B_kj=&(T->val[kj_ptrB*B_block_size]);
434             for (irb=0; irb<row_block_size; ++irb) {                     for (irb=0; irb<row_block_size; ++irb) {
435               for (icb=0; icb<col_block_size; ++icb) {                       for (icb=0; icb<col_block_size; ++icb) {
436              rtmp=C_ij[irb+row_block_size*icb];                          rtmp=C_ij[irb+row_block_size*icb];
437              for (ib=0; ib<A_col_block_size; ++ib) {                          for (ib=0; ib<A_col_block_size; ++ib) {
438                rtmp+=A_ik[irb+row_block_size*ib]*B_kj[ib+A_col_block_size*icb];                            rtmp+=A_ik[irb+row_block_size*ib]*B_kj[ib+A_col_block_size*icb];
439              }                          }
440              C_ij[irb+row_block_size*icb]=rtmp;                          C_ij[irb+row_block_size*icb]=rtmp;
441               }                       }
442                     }                     }
443                     ik_ptrA ++;                     ik_ptrA ++;
444                     kj_ptrB ++;                     kj_ptrB ++;
# Line 440  void SparseMatrix_MatrixMatrixTranspose_ Line 447  void SparseMatrix_MatrixMatrixTranspose_
447                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
448                   }                   }
449                 }                 }
450          }              }
451       }           }
452        } /* end of parallel region */        } /* end of parallel region */
453                
454     }     }
455  }  }
456    
457  /* not good for block size 1 */  /* not good for block size 1 */
458  void SparseMatrix_MatrixMatrixTranspose_DB(SparseMatrix* C, const SparseMatrix* A, const SparseMatrix* B, const SparseMatrix* T)  void SparseMatrix_MatrixMatrixTranspose_DB(SparseMatrix_ptr C, const_SparseMatrix_ptr A, const_SparseMatrix_ptr B, const_SparseMatrix_ptr T)
459  {  {
460     const dim_t n = C->numRows;     const dim_t n = C->numRows;
461     const dim_t row_block_size = C->row_block_size;     const dim_t row_block_size = C->row_block_size;
# Line 465  void SparseMatrix_MatrixMatrixTranspose_ Line 472  void SparseMatrix_MatrixMatrixTranspose_
472     if ( (row_block_size == 2) && (col_block_size ==2 ) && (A_block_size == 2) ) {     if ( (row_block_size == 2) && (col_block_size ==2 ) && (A_block_size == 2) ) {
473        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_01, C_ij_11)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_01, C_ij_11)
474        {        {
475       #pragma omp for schedule(static)           #pragma omp for schedule(static)
476       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
477          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
478             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
479                            
480             C_ij_00=0;                 C_ij_00=0;
481             C_ij_10=0;                 C_ij_10=0;
482             C_ij_01=0;                 C_ij_01=0;
483             C_ij_11=0;                 C_ij_11=0;
484    
485             C_ij=&(C->val[ij_ptrC*4]);                 C_ij=&(C->val[ij_ptrC*4]);
486    
487             ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
488                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
489                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
490                 kjb = T->pattern->ptr[j+1];                 kjb = T->pattern->ptr[j+1];
# Line 493  void SparseMatrix_MatrixMatrixTranspose_ Line 500  void SparseMatrix_MatrixMatrixTranspose_
500                     if (kj_ptrB >= kjb) break;                     if (kj_ptrB >= kjb) break;
501                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
502                   } else {                   } else {
503              A_ik=&(A->val[ik_ptrA*2]);                          A_ik=&(A->val[ik_ptrA*2]);
504              B_kj=&(T->val[kj_ptrB*4]);                          B_kj=&(T->val[kj_ptrB*4]);
505                                        
506              C_ij_00 +=A_ik[0]*B_kj[0+2*0];                          C_ij_00 +=A_ik[0]*B_kj[0+2*0];
507              C_ij_10 +=A_ik[1]*B_kj[1+2*0];                          C_ij_10 +=A_ik[1]*B_kj[1+2*0];
508                                                    
509              C_ij_01 +=A_ik[0]*B_kj[0+2*1];                          C_ij_01 +=A_ik[0]*B_kj[0+2*1];
510              C_ij_11 +=A_ik[1]*B_kj[1+2*1];                          C_ij_11 +=A_ik[1]*B_kj[1+2*1];
511              ik_ptrA ++;                          ik_ptrA ++;
512              kj_ptrB ++;                          kj_ptrB ++;
513              if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                          if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
514              kA=A->pattern->index[ik_ptrA];                          kA=A->pattern->index[ik_ptrA];
515              kB=T->pattern->index[kj_ptrB];                          kB=T->pattern->index[kj_ptrB];
   
          }  
            }  
            C_ij[0+2*0]=C_ij_00;  
            C_ij[1+2*0]=C_ij_10;  
            C_ij[0+2*1]=C_ij_01;  
            C_ij[1+2*1]=C_ij_11;  
516    
517          }                   }
518       }                 }
519                   C_ij[0+2*0]=C_ij_00;
520                   C_ij[1+2*0]=C_ij_10;
521                   C_ij[0+2*1]=C_ij_01;
522                   C_ij[1+2*1]=C_ij_11;
523    
524                }
525             }
526        } /* end of parallel region */        } /* end of parallel region */
527                
528     } else if ( (row_block_size == 3) && (col_block_size ==3 ) && (A_block_size == 3)  ){     } else if ( (row_block_size == 3) && (col_block_size ==3 ) && (A_block_size == 3)  ){
529        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_01, C_ij_11, C_ij_21, C_ij_02, C_ij_12, C_ij_22)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_01, C_ij_11, C_ij_21, C_ij_02, C_ij_12, C_ij_22)
530        {        {
531       #pragma omp for schedule(static)           #pragma omp for schedule(static)
532       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
533          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
534             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
535                            
536             C_ij_00=0;                 C_ij_00=0;
537             C_ij_10=0;                 C_ij_10=0;
538             C_ij_20=0;                 C_ij_20=0;
539             C_ij_01=0;                 C_ij_01=0;
540             C_ij_11=0;                 C_ij_11=0;
541             C_ij_21=0;                 C_ij_21=0;
542             C_ij_02=0;                 C_ij_02=0;
543             C_ij_12=0;                 C_ij_12=0;
544             C_ij_22=0;                 C_ij_22=0;
545                            
546             C_ij=&(C->val[ij_ptrC*9]);                 C_ij=&(C->val[ij_ptrC*9]);
547    
548             ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
549                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
550                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
551                 kjb = T->pattern->ptr[j+1];                 kjb = T->pattern->ptr[j+1];
# Line 554  void SparseMatrix_MatrixMatrixTranspose_ Line 561  void SparseMatrix_MatrixMatrixTranspose_
561                     if (kj_ptrB >= kjb) break;                     if (kj_ptrB >= kjb) break;
562                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
563                   } else {                   } else {
564              A_ik=&(A->val[ik_ptrA*3]);                          A_ik=&(A->val[ik_ptrA*3]);
565              B_kj=&(T->val[kj_ptrB*9]);                          B_kj=&(T->val[kj_ptrB*9]);
566                                        
567              C_ij_00 +=A_ik[0]*B_kj[0+3*0];                          C_ij_00 +=A_ik[0]*B_kj[0+3*0];
568              C_ij_10 +=A_ik[1]*B_kj[1+3*0];                          C_ij_10 +=A_ik[1]*B_kj[1+3*0];
569              C_ij_20 +=A_ik[2]*B_kj[2+3*0];                          C_ij_20 +=A_ik[2]*B_kj[2+3*0];
570                                        
571              C_ij_01 +=A_ik[0]*B_kj[0+3*1];                          C_ij_01 +=A_ik[0]*B_kj[0+3*1];
572              C_ij_11 +=A_ik[1]*B_kj[1+3*1];                          C_ij_11 +=A_ik[1]*B_kj[1+3*1];
573              C_ij_21 +=A_ik[2]*B_kj[2+3*1];                          C_ij_21 +=A_ik[2]*B_kj[2+3*1];
574                                        
575              C_ij_02 +=A_ik[0]*B_kj[0+3*2];                          C_ij_02 +=A_ik[0]*B_kj[0+3*2];
576              C_ij_12 +=A_ik[1]*B_kj[1+3*2];                          C_ij_12 +=A_ik[1]*B_kj[1+3*2];
577              C_ij_22 +=A_ik[2]*B_kj[2+3*2];                          C_ij_22 +=A_ik[2]*B_kj[2+3*2];
578              ik_ptrA ++;                          ik_ptrA ++;
579              kj_ptrB ++;                          kj_ptrB ++;
580              if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                          if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
581              kA=A->pattern->index[ik_ptrA];                          kA=A->pattern->index[ik_ptrA];
582              kB=T->pattern->index[kj_ptrB];                          kB=T->pattern->index[kj_ptrB];
583           }                   }
584             }                 }
585             C_ij[0+3*0]=C_ij_00;                 C_ij[0+3*0]=C_ij_00;
586             C_ij[1+3*0]=C_ij_10;                 C_ij[1+3*0]=C_ij_10;
587             C_ij[2+3*0]=C_ij_20;                 C_ij[2+3*0]=C_ij_20;
588             C_ij[0+3*1]=C_ij_01;                 C_ij[0+3*1]=C_ij_01;
589             C_ij[1+3*1]=C_ij_11;                 C_ij[1+3*1]=C_ij_11;
590             C_ij[2+3*1]=C_ij_21;                 C_ij[2+3*1]=C_ij_21;
591             C_ij[0+3*2]=C_ij_02;                 C_ij[0+3*2]=C_ij_02;
592             C_ij[1+3*2]=C_ij_12;                 C_ij[1+3*2]=C_ij_12;
593             C_ij[2+3*2]=C_ij_22;                 C_ij[2+3*2]=C_ij_22;
594          }              }
595       }           }
596        } /* end of parallel region */        } /* end of parallel region */
597     } else if ( (row_block_size == 4) && (col_block_size ==4 ) && (A_block_size == 4)  ){     } else if ( (row_block_size == 4) && (col_block_size ==4 ) && (A_block_size == 4)  ){
598        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_30, C_ij_01, C_ij_11, C_ij_21, C_ij_31, C_ij_02, C_ij_12, C_ij_22, C_ij_32, C_ij_03, C_ij_13, C_ij_23, C_ij_33)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_30, C_ij_01, C_ij_11, C_ij_21, C_ij_31, C_ij_02, C_ij_12, C_ij_22, C_ij_32, C_ij_03, C_ij_13, C_ij_23, C_ij_33)
599        {        {
600       #pragma omp for schedule(static)           #pragma omp for schedule(static)
601       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
602          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
603             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
604                            
605             C_ij_00=0;                 C_ij_00=0;
606             C_ij_10=0;                 C_ij_10=0;
607             C_ij_20=0;                 C_ij_20=0;
608             C_ij_30=0;                 C_ij_30=0;
609             C_ij_01=0;                 C_ij_01=0;
610             C_ij_11=0;                 C_ij_11=0;
611             C_ij_21=0;                 C_ij_21=0;
612             C_ij_31=0;                 C_ij_31=0;
613             C_ij_02=0;                 C_ij_02=0;
614             C_ij_12=0;                 C_ij_12=0;
615             C_ij_22=0;                 C_ij_22=0;
616             C_ij_32=0;                 C_ij_32=0;
617             C_ij_03=0;                 C_ij_03=0;
618             C_ij_13=0;                 C_ij_13=0;
619             C_ij_23=0;                 C_ij_23=0;
620             C_ij_33=0;                           C_ij_33=0;              
621                            
622             C_ij=&(C->val[ij_ptrC*16]);                 C_ij=&(C->val[ij_ptrC*16]);
623    
624             ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
625                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
626                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
627                 kjb = T->pattern->ptr[j+1];                 kjb = T->pattern->ptr[j+1];
# Line 630  void SparseMatrix_MatrixMatrixTranspose_ Line 637  void SparseMatrix_MatrixMatrixTranspose_
637                     if (kj_ptrB >= kjb) break;                     if (kj_ptrB >= kjb) break;
638                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
639                   } else {                   } else {
640                 A_ik=&(A->val[ik_ptrA*4]);                             A_ik=&(A->val[ik_ptrA*4]);
641                 B_kj=&(T->val[kj_ptrB*16]);                             B_kj=&(T->val[kj_ptrB*16]);
642                                            
643                 C_ij_00 +=A_ik[0]*B_kj[0+4*0];                             C_ij_00 +=A_ik[0]*B_kj[0+4*0];
644                 C_ij_10 +=A_ik[1]*B_kj[1+4*0];                             C_ij_10 +=A_ik[1]*B_kj[1+4*0];
645                 C_ij_20 +=A_ik[2]*B_kj[2+4*0];                             C_ij_20 +=A_ik[2]*B_kj[2+4*0];
646                 C_ij_30 +=A_ik[3]*B_kj[3+4*0];                             C_ij_30 +=A_ik[3]*B_kj[3+4*0];
647                                                            
648                 C_ij_01 +=A_ik[0]*B_kj[0+4*1];                             C_ij_01 +=A_ik[0]*B_kj[0+4*1];
649                 C_ij_11 +=A_ik[1]*B_kj[1+4*1];                             C_ij_11 +=A_ik[1]*B_kj[1+4*1];
650                 C_ij_21 +=A_ik[2]*B_kj[2+4*1];                             C_ij_21 +=A_ik[2]*B_kj[2+4*1];
651                 C_ij_31 +=A_ik[3]*B_kj[3+4*1];                             C_ij_31 +=A_ik[3]*B_kj[3+4*1];
652                                            
653                 C_ij_02 +=A_ik[0]*B_kj[0+4*2];                             C_ij_02 +=A_ik[0]*B_kj[0+4*2];
654                 C_ij_12 +=A_ik[1]*B_kj[1+4*2];                             C_ij_12 +=A_ik[1]*B_kj[1+4*2];
655                 C_ij_22 +=A_ik[2]*B_kj[2+4*2];                             C_ij_22 +=A_ik[2]*B_kj[2+4*2];
656                 C_ij_32 +=A_ik[3]*B_kj[3+4*2];                             C_ij_32 +=A_ik[3]*B_kj[3+4*2];
657                                            
658                 C_ij_03 +=A_ik[0]*B_kj[0+4*3];                             C_ij_03 +=A_ik[0]*B_kj[0+4*3];
659                 C_ij_13 +=A_ik[1]*B_kj[1+4*3];                             C_ij_13 +=A_ik[1]*B_kj[1+4*3];
660                 C_ij_23 +=A_ik[2]*B_kj[2+4*3];                             C_ij_23 +=A_ik[2]*B_kj[2+4*3];
661                 C_ij_33 +=A_ik[3]*B_kj[3+4*3];                             C_ij_33 +=A_ik[3]*B_kj[3+4*3];
662    
663                 ik_ptrA ++;                             ik_ptrA ++;
664                 kj_ptrB ++;                             kj_ptrB ++;
665                 if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                             if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
666                 kA=A->pattern->index[ik_ptrA];                             kA=A->pattern->index[ik_ptrA];
667                 kB=T->pattern->index[kj_ptrB];                             kB=T->pattern->index[kj_ptrB];
668           }                   }
669             }                 }
670             C_ij[0+4*0]=C_ij_00;                 C_ij[0+4*0]=C_ij_00;
671             C_ij[1+4*0]=C_ij_10;                 C_ij[1+4*0]=C_ij_10;
672             C_ij[2+4*0]=C_ij_20;                 C_ij[2+4*0]=C_ij_20;
673             C_ij[3+4*0]=C_ij_30;                 C_ij[3+4*0]=C_ij_30;
674             C_ij[0+4*1]=C_ij_01;                 C_ij[0+4*1]=C_ij_01;
675             C_ij[1+4*1]=C_ij_11;                 C_ij[1+4*1]=C_ij_11;
676             C_ij[2+4*1]=C_ij_21;                 C_ij[2+4*1]=C_ij_21;
677             C_ij[3+4*1]=C_ij_31;                 C_ij[3+4*1]=C_ij_31;
678             C_ij[0+4*2]=C_ij_02;                 C_ij[0+4*2]=C_ij_02;
679             C_ij[1+4*2]=C_ij_12;                 C_ij[1+4*2]=C_ij_12;
680             C_ij[2+4*2]=C_ij_22;                 C_ij[2+4*2]=C_ij_22;
681             C_ij[3+4*2]=C_ij_32;                 C_ij[3+4*2]=C_ij_32;    
682             C_ij[0+4*3]=C_ij_03;                 C_ij[0+4*3]=C_ij_03;
683             C_ij[1+4*3]=C_ij_13;                 C_ij[1+4*3]=C_ij_13;
684             C_ij[2+4*3]=C_ij_23;                 C_ij[2+4*3]=C_ij_23;
685             C_ij[3+4*3]=C_ij_33;                       C_ij[3+4*3]=C_ij_33;              
686          }              }
687       }           }
688        } /* end of parallel region */        } /* end of parallel region */
689                    
690     } else {     } else {
691        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, rtmp, A_ik,B_kj )        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, rtmp, A_ik,B_kj )
692        {        {
693       #pragma omp for schedule(static)           #pragma omp for schedule(static)
694       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
695          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
696             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
697             C_ij=&(C->val[ij_ptrC*C_block_size]);                 C_ij=&(C->val[ij_ptrC*C_block_size]);
698             for (ib=0; ib<C_block_size; ++ib)  C_ij[ib]=0;                 for (ib=0; ib<C_block_size; ++ib)  C_ij[ib]=0;
699             ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
700                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
701                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
702                 kjb = T->pattern->ptr[j+1];                 kjb = T->pattern->ptr[j+1];
# Line 705  void SparseMatrix_MatrixMatrixTranspose_ Line 712  void SparseMatrix_MatrixMatrixTranspose_
712                     if (kj_ptrB >= kjb) break;                     if (kj_ptrB >= kjb) break;
713                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
714                   } else {                   } else {
715                 A_ik=&(A->val[ik_ptrA*A_block_size]);                         A_ik=&(A->val[ik_ptrA*A_block_size]);
716                 B_kj=&(T->val[kj_ptrB*B_block_size]);                         B_kj=&(T->val[kj_ptrB*B_block_size]);
717                                        
718                 for (irb=0; irb<A_block_size; ++irb) {                         for (irb=0; irb<A_block_size; ++irb) {
719                rtmp=A_ik[irb];                            rtmp=A_ik[irb];
720                for (icb=0; icb<col_block_size; ++icb) {                            for (icb=0; icb<col_block_size; ++icb) {
721                      C_ij[irb+row_block_size*icb]+=rtmp*B_kj[irb+A_col_block_size*icb];                                  C_ij[irb+row_block_size*icb]+=rtmp*B_kj[irb+A_col_block_size*icb];
722                }                            }
723                 }                         }
724                 ik_ptrA ++;                         ik_ptrA ++;
725                 kj_ptrB ++;                         kj_ptrB ++;
726                 if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                         if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
727                 kA=A->pattern->index[ik_ptrA];                         kA=A->pattern->index[ik_ptrA];
728                 kB=T->pattern->index[kj_ptrB];                         kB=T->pattern->index[kj_ptrB];
729            }                    }
730             }                 }
731          }              }
732       }           }
733        } /* end of parallel region */        } /* end of parallel region */
734                
735     }     }
736  }  }
737    
738  /* not good for block size 1 */  /* not good for block size 1 */
739  void SparseMatrix_MatrixMatrixTranspose_BD(SparseMatrix* C, const SparseMatrix* A, const SparseMatrix* B, const SparseMatrix* T)  void SparseMatrix_MatrixMatrixTranspose_BD(SparseMatrix_ptr C, const_SparseMatrix_ptr A, const_SparseMatrix_ptr B, const_SparseMatrix_ptr T)
740  {  {
741     const dim_t n = C->numRows;     const dim_t n = C->numRows;
742     const dim_t row_block_size = C->row_block_size;     const dim_t row_block_size = C->row_block_size;
# Line 745  void SparseMatrix_MatrixMatrixTranspose_ Line 752  void SparseMatrix_MatrixMatrixTranspose_
752     if ( (row_block_size == 2) && (col_block_size ==2 ) && (B_block_size == 2) ) {     if ( (row_block_size == 2) && (col_block_size ==2 ) && (B_block_size == 2) ) {
753        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_01, C_ij_11)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_01, C_ij_11)
754        {        {
755       #pragma omp for schedule(static)           #pragma omp for schedule(static)
756       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
757          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
758             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
759                            
760             C_ij_00=0;                 C_ij_00=0;
761             C_ij_10=0;                 C_ij_10=0;
762             C_ij_01=0;                 C_ij_01=0;
763             C_ij_11=0;                 C_ij_11=0;
764    
765             C_ij=&(C->val[ij_ptrC*4]);                 C_ij=&(C->val[ij_ptrC*4]);
766    
767             ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
768                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
769                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
770                 kjb = T->pattern->ptr[j+1];                 kjb = T->pattern->ptr[j+1];
# Line 773  void SparseMatrix_MatrixMatrixTranspose_ Line 780  void SparseMatrix_MatrixMatrixTranspose_
780                     if (kj_ptrB >= kjb) break;                     if (kj_ptrB >= kjb) break;
781                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
782                   } else {                   } else {
783              A_ik=&(A->val[ik_ptrA*4]);                          A_ik=&(A->val[ik_ptrA*4]);
784              B_kj=&(T->val[kj_ptrB*2]);                          B_kj=&(T->val[kj_ptrB*2]);
785                                        
786              C_ij_00 +=A_ik[0+2*0]*B_kj[0];                          C_ij_00 +=A_ik[0+2*0]*B_kj[0];
787              C_ij_10 +=A_ik[1+2*0]*B_kj[0];                          C_ij_10 +=A_ik[1+2*0]*B_kj[0];
788                                                    
789              C_ij_01 +=A_ik[0+2*1]*B_kj[1];                          C_ij_01 +=A_ik[0+2*1]*B_kj[1];
790              C_ij_11 +=A_ik[1+2*1]*B_kj[1];                          C_ij_11 +=A_ik[1+2*1]*B_kj[1];
791    
792              ik_ptrA ++;                          ik_ptrA ++;
793              kj_ptrB ++;                          kj_ptrB ++;
794              if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                          if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
795              kA=A->pattern->index[ik_ptrA];                          kA=A->pattern->index[ik_ptrA];
796              kB=T->pattern->index[kj_ptrB];                          kB=T->pattern->index[kj_ptrB];
797           }                   }
798             }                 }
799             C_ij[0+2*0]=C_ij_00;                 C_ij[0+2*0]=C_ij_00;
800             C_ij[1+2*0]=C_ij_10;                 C_ij[1+2*0]=C_ij_10;
801             C_ij[0+2*1]=C_ij_01;                 C_ij[0+2*1]=C_ij_01;
802             C_ij[1+2*1]=C_ij_11;                 C_ij[1+2*1]=C_ij_11;
803    
804          }              }
805       }           }
806        } /* end of parallel region */        } /* end of parallel region */
807                
808     } else if ( (row_block_size == 3) && (col_block_size ==3 ) && (B_block_size == 3)  ){     } else if ( (row_block_size == 3) && (col_block_size ==3 ) && (B_block_size == 3)  ){
809        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_01, C_ij_11, C_ij_21, C_ij_02, C_ij_12, C_ij_22)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_01, C_ij_11, C_ij_21, C_ij_02, C_ij_12, C_ij_22)
810        {        {
811       #pragma omp for schedule(static)           #pragma omp for schedule(static)
812       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
813          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
814             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
815                            
816             C_ij_00=0;                 C_ij_00=0;
817             C_ij_10=0;                 C_ij_10=0;
818             C_ij_20=0;                 C_ij_20=0;
819             C_ij_01=0;                 C_ij_01=0;
820             C_ij_11=0;                 C_ij_11=0;
821             C_ij_21=0;                 C_ij_21=0;
822             C_ij_02=0;                 C_ij_02=0;
823             C_ij_12=0;                 C_ij_12=0;
824             C_ij_22=0;                 C_ij_22=0;
825                            
826             C_ij=&(C->val[ij_ptrC*9]);                 C_ij=&(C->val[ij_ptrC*9]);
827    
828             ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
829                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
830                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
831                 kjb = T->pattern->ptr[j+1];                 kjb = T->pattern->ptr[j+1];
# Line 834  void SparseMatrix_MatrixMatrixTranspose_ Line 841  void SparseMatrix_MatrixMatrixTranspose_
841                     if (kj_ptrB >= kjb) break;                     if (kj_ptrB >= kjb) break;
842                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
843                   } else {                   } else {
844              A_ik=&(A->val[ik_ptrA*9]);                          A_ik=&(A->val[ik_ptrA*9]);
845              B_kj=&(T->val[kj_ptrB*3]);                          B_kj=&(T->val[kj_ptrB*3]);
846                                        
847              C_ij_00 +=A_ik[0+3*0]*B_kj[0];                          C_ij_00 +=A_ik[0+3*0]*B_kj[0];
848              C_ij_10 +=A_ik[1+3*0]*B_kj[0];                          C_ij_10 +=A_ik[1+3*0]*B_kj[0];
849              C_ij_20 +=A_ik[2+3*0]*B_kj[0];                          C_ij_20 +=A_ik[2+3*0]*B_kj[0];
850                                        
851              C_ij_01 +=A_ik[0+3*1]*B_kj[1];                          C_ij_01 +=A_ik[0+3*1]*B_kj[1];
852              C_ij_11 +=A_ik[1+3*1]*B_kj[1];                          C_ij_11 +=A_ik[1+3*1]*B_kj[1];
853              C_ij_21 +=A_ik[2+3*1]*B_kj[1];                          C_ij_21 +=A_ik[2+3*1]*B_kj[1];
854    
855              C_ij_02 +=A_ik[0+3*2]*B_kj[2];                          C_ij_02 +=A_ik[0+3*2]*B_kj[2];
856              C_ij_12 +=A_ik[1+3*2]*B_kj[2];                          C_ij_12 +=A_ik[1+3*2]*B_kj[2];
857              C_ij_22 +=A_ik[2+3*2]*B_kj[2];                          C_ij_22 +=A_ik[2+3*2]*B_kj[2];
858    
859              ik_ptrA ++;                          ik_ptrA ++;
860              kj_ptrB ++;                          kj_ptrB ++;
861              if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                          if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
862              kA=A->pattern->index[ik_ptrA];                          kA=A->pattern->index[ik_ptrA];
863              kB=T->pattern->index[kj_ptrB];                          kB=T->pattern->index[kj_ptrB];
864           }                   }
865             }                 }
866             C_ij[0+3*0]=C_ij_00;                 C_ij[0+3*0]=C_ij_00;
867             C_ij[1+3*0]=C_ij_10;                 C_ij[1+3*0]=C_ij_10;
868             C_ij[2+3*0]=C_ij_20;                 C_ij[2+3*0]=C_ij_20;
869             C_ij[0+3*1]=C_ij_01;                 C_ij[0+3*1]=C_ij_01;
870             C_ij[1+3*1]=C_ij_11;                 C_ij[1+3*1]=C_ij_11;
871             C_ij[2+3*1]=C_ij_21;                 C_ij[2+3*1]=C_ij_21;
872             C_ij[0+3*2]=C_ij_02;                 C_ij[0+3*2]=C_ij_02;
873             C_ij[1+3*2]=C_ij_12;                 C_ij[1+3*2]=C_ij_12;
874             C_ij[2+3*2]=C_ij_22;                 C_ij[2+3*2]=C_ij_22;
875          }              }
876       }           }
877        } /* end of parallel region */        } /* end of parallel region */
878     } else if ( (row_block_size == 4) && (col_block_size ==4 ) && (B_block_size == 4)  ){     } else if ( (row_block_size == 4) && (col_block_size ==4 ) && (B_block_size == 4)  ){
879        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_30, C_ij_01, C_ij_11, C_ij_21, C_ij_31, C_ij_02, C_ij_12, C_ij_22, C_ij_32, C_ij_03, C_ij_13, C_ij_23, C_ij_33)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, A_ik,B_kj,C_ij_00, C_ij_10, C_ij_20, C_ij_30, C_ij_01, C_ij_11, C_ij_21, C_ij_31, C_ij_02, C_ij_12, C_ij_22, C_ij_32, C_ij_03, C_ij_13, C_ij_23, C_ij_33)
880        {        {
881       #pragma omp for schedule(static)           #pragma omp for schedule(static)
882       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
883          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
884             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
885                            
886             C_ij_00=0;                 C_ij_00=0;
887             C_ij_10=0;                 C_ij_10=0;
888             C_ij_20=0;                 C_ij_20=0;
889             C_ij_30=0;                 C_ij_30=0;
890             C_ij_01=0;                 C_ij_01=0;
891             C_ij_11=0;                 C_ij_11=0;
892             C_ij_21=0;                 C_ij_21=0;
893             C_ij_31=0;                 C_ij_31=0;
894             C_ij_02=0;                 C_ij_02=0;
895             C_ij_12=0;                 C_ij_12=0;
896             C_ij_22=0;                 C_ij_22=0;
897             C_ij_32=0;                 C_ij_32=0;
898             C_ij_03=0;                 C_ij_03=0;
899             C_ij_13=0;                 C_ij_13=0;
900             C_ij_23=0;                 C_ij_23=0;
901             C_ij_33=0;                           C_ij_33=0;              
902                            
903             C_ij=&(C->val[ij_ptrC*16]);                 C_ij=&(C->val[ij_ptrC*16]);
904    
905             ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
906                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
907                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
908                 kjb = T->pattern->ptr[j+1];                 kjb = T->pattern->ptr[j+1];
# Line 911  void SparseMatrix_MatrixMatrixTranspose_ Line 918  void SparseMatrix_MatrixMatrixTranspose_
918                     if (kj_ptrB >= kjb) break;                     if (kj_ptrB >= kjb) break;
919                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
920                   } else {                   } else {
921                 A_ik=&(A->val[ik_ptrA*16]);                             A_ik=&(A->val[ik_ptrA*16]);
922                 B_kj=&(T->val[kj_ptrB*4]);                             B_kj=&(T->val[kj_ptrB*4]);
923                                            
924                 C_ij_00 +=A_ik[0+4*0]*B_kj[0];                             C_ij_00 +=A_ik[0+4*0]*B_kj[0];
925                 C_ij_10 +=A_ik[1+4*0]*B_kj[0];                             C_ij_10 +=A_ik[1+4*0]*B_kj[0];
926                 C_ij_20 +=A_ik[2+4*0]*B_kj[0];                             C_ij_20 +=A_ik[2+4*0]*B_kj[0];
927                 C_ij_30 +=A_ik[3+4*0]*B_kj[0];                             C_ij_30 +=A_ik[3+4*0]*B_kj[0];
928                                                            
929                 C_ij_01 +=A_ik[0+4*1]*B_kj[1];                             C_ij_01 +=A_ik[0+4*1]*B_kj[1];
930                 C_ij_11 +=A_ik[1+4*1]*B_kj[1];                             C_ij_11 +=A_ik[1+4*1]*B_kj[1];
931                 C_ij_21 +=A_ik[2+4*1]*B_kj[1];                             C_ij_21 +=A_ik[2+4*1]*B_kj[1];
932                 C_ij_31 +=A_ik[3+4*1]*B_kj[1];                             C_ij_31 +=A_ik[3+4*1]*B_kj[1];
933    
934                 C_ij_02 +=A_ik[0+4*2]*B_kj[2];                             C_ij_02 +=A_ik[0+4*2]*B_kj[2];
935                 C_ij_12 +=A_ik[1+4*2]*B_kj[2];                             C_ij_12 +=A_ik[1+4*2]*B_kj[2];
936                 C_ij_22 +=A_ik[2+4*2]*B_kj[2];                             C_ij_22 +=A_ik[2+4*2]*B_kj[2];
937                 C_ij_32 +=A_ik[3+4*2]*B_kj[2];                             C_ij_32 +=A_ik[3+4*2]*B_kj[2];
938    
939                 C_ij_03 +=A_ik[0+4*3]*B_kj[3];                             C_ij_03 +=A_ik[0+4*3]*B_kj[3];
940                 C_ij_13 +=A_ik[1+4*3]*B_kj[3];                             C_ij_13 +=A_ik[1+4*3]*B_kj[3];
941                 C_ij_23 +=A_ik[2+4*3]*B_kj[3];                             C_ij_23 +=A_ik[2+4*3]*B_kj[3];
942                 C_ij_33 +=A_ik[3+4*3]*B_kj[3];                             C_ij_33 +=A_ik[3+4*3]*B_kj[3];
943    
944                 ik_ptrA ++;                             ik_ptrA ++;
945                 kj_ptrB ++;                             kj_ptrB ++;
946                 if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                             if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
947                 kA=A->pattern->index[ik_ptrA];                             kA=A->pattern->index[ik_ptrA];
948                 kB=T->pattern->index[kj_ptrB];                             kB=T->pattern->index[kj_ptrB];
949           }                   }
950             }                 }
951             C_ij[0+4*0]=C_ij_00;                 C_ij[0+4*0]=C_ij_00;
952             C_ij[1+4*0]=C_ij_10;                 C_ij[1+4*0]=C_ij_10;
953             C_ij[2+4*0]=C_ij_20;                 C_ij[2+4*0]=C_ij_20;
954             C_ij[3+4*0]=C_ij_30;                 C_ij[3+4*0]=C_ij_30;
955             C_ij[0+4*1]=C_ij_01;                 C_ij[0+4*1]=C_ij_01;
956             C_ij[1+4*1]=C_ij_11;                 C_ij[1+4*1]=C_ij_11;
957             C_ij[2+4*1]=C_ij_21;                 C_ij[2+4*1]=C_ij_21;
958             C_ij[3+4*1]=C_ij_31;                 C_ij[3+4*1]=C_ij_31;
959             C_ij[0+4*2]=C_ij_02;                 C_ij[0+4*2]=C_ij_02;
960             C_ij[1+4*2]=C_ij_12;                 C_ij[1+4*2]=C_ij_12;
961             C_ij[2+4*2]=C_ij_22;                 C_ij[2+4*2]=C_ij_22;
962             C_ij[3+4*2]=C_ij_32;                 C_ij[3+4*2]=C_ij_32;
963             C_ij[0+4*3]=C_ij_03;                 C_ij[0+4*3]=C_ij_03;
964             C_ij[1+4*3]=C_ij_13;                 C_ij[1+4*3]=C_ij_13;
965             C_ij[2+4*3]=C_ij_23;                 C_ij[2+4*3]=C_ij_23;
966             C_ij[3+4*3]=C_ij_33;                       C_ij[3+4*3]=C_ij_33;              
967          }              }
968       }           }
969        } /* end of parallel region */        } /* end of parallel region */
970                    
971     } else {     } else {
972        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, rtmp, A_ik,B_kj )        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, irb, icb, ib, rtmp, A_ik,B_kj )
973        {        {
974       #pragma omp for schedule(static)           #pragma omp for schedule(static)
975       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
976          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
977             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
978             C_ij=&(C->val[ij_ptrC*C_block_size]);                 C_ij=&(C->val[ij_ptrC*C_block_size]);
979             for (ib=0; ib<C_block_size; ++ib)  C_ij[ib]=0;                 for (ib=0; ib<C_block_size; ++ib)  C_ij[ib]=0;
980             ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
981                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
982                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
983                 kjb = T->pattern->ptr[j+1];                 kjb = T->pattern->ptr[j+1];
# Line 986  void SparseMatrix_MatrixMatrixTranspose_ Line 993  void SparseMatrix_MatrixMatrixTranspose_
993                     if (kj_ptrB >= kjb) break;                     if (kj_ptrB >= kjb) break;
994                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
995                   } else {                   } else {
996                 A_ik=&(A->val[ik_ptrA*A_block_size]);                         A_ik=&(A->val[ik_ptrA*A_block_size]);
997                 B_kj=&(T->val[kj_ptrB*B_block_size]);                         B_kj=&(T->val[kj_ptrB*B_block_size]);
998    
999                 for (icb=0; icb<B_block_size; ++icb) {                         for (icb=0; icb<B_block_size; ++icb) {
1000                rtmp=B_kj[icb];                            rtmp=B_kj[icb];
1001                for (irb=0; irb<row_block_size; ++irb) {                            for (irb=0; irb<row_block_size; ++irb) {
1002                   C_ij[irb+row_block_size*icb]+=A_ik[irb+row_block_size*icb]*rtmp;                               C_ij[irb+row_block_size*icb]+=A_ik[irb+row_block_size*icb]*rtmp;
1003                }                            }
1004                 }                         }
1005    
1006                 ik_ptrA ++;                         ik_ptrA ++;
1007                 kj_ptrB ++;                         kj_ptrB ++;
1008                 if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                         if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
1009                 kA=A->pattern->index[ik_ptrA];                         kA=A->pattern->index[ik_ptrA];
1010                 kB=T->pattern->index[kj_ptrB];                         kB=T->pattern->index[kj_ptrB];
1011           }                   }
1012             }                 }
1013          }              }
1014       }           }
1015        } /* end of parallel region */        } /* end of parallel region */
1016                
1017     }     }
1018  }  }
1019    
1020  /* not good for block size 1 */  /* not good for block size 1 */
1021  void SparseMatrix_MatrixMatrixTranspose_DD(SparseMatrix* C, const SparseMatrix* A, const SparseMatrix* B, const SparseMatrix* T)  void SparseMatrix_MatrixMatrixTranspose_DD(SparseMatrix_ptr C, const_SparseMatrix_ptr A, const_SparseMatrix_ptr B, const_SparseMatrix_ptr T)
1022  {  {
1023     const dim_t n = C->numRows;     const dim_t n = C->numRows;
1024     const dim_t C_block_size =C->block_size;     const dim_t C_block_size =C->block_size;
# Line 1025  void SparseMatrix_MatrixMatrixTranspose_ Line 1032  void SparseMatrix_MatrixMatrixTranspose_
1032     if ( (A_block_size == 1) && (B_block_size ==1) && (C_block_size == 1) ) {     if ( (A_block_size == 1) && (B_block_size ==1) && (C_block_size == 1) ) {
1033        #pragma omp parallel private(i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, C_ij_0)        #pragma omp parallel private(i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, C_ij_0)
1034        {        {
1035       #pragma omp for schedule(static)           #pragma omp for schedule(static)
1036           for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
1037          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
1038             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
1039             C_ij_0=0;                 C_ij_0=0;
1040                 ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
1041                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
1042                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
# Line 1054  void SparseMatrix_MatrixMatrixTranspose_ Line 1061  void SparseMatrix_MatrixMatrixTranspose_
1061                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
1062                   }                   }
1063                 }                 }
1064             C->val[ij_ptrC]=C_ij_0;                 C->val[ij_ptrC]=C_ij_0;
1065          }              }
1066           }           }
1067        } /* end of parallel region */        } /* end of parallel region */
1068      } else if ( (A_block_size == 2) && (B_block_size ==2) && (C_block_size == 2) ) {      } else if ( (A_block_size == 2) && (B_block_size ==2) && (C_block_size == 2) ) {
1069        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, A_ik,B_kj, C_ij_0, C_ij_1)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, A_ik,B_kj, C_ij_0, C_ij_1)
1070        {        {
1071       #pragma omp for schedule(static)           #pragma omp for schedule(static)
1072       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
1073          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
1074             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
1075             C_ij=&(C->val[ij_ptrC*2]);                 C_ij=&(C->val[ij_ptrC*2]);
1076             C_ij_0=0;                 C_ij_0=0;
1077             C_ij_1=0;                 C_ij_1=0;
1078                 ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
1079                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
1080                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
# Line 1095  void SparseMatrix_MatrixMatrixTranspose_ Line 1102  void SparseMatrix_MatrixMatrixTranspose_
1102                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
1103                   }                   }
1104                 }                 }
1105             C_ij[0]=C_ij_0;                 C_ij[0]=C_ij_0;
1106             C_ij[1]=C_ij_1;                 C_ij[1]=C_ij_1;
1107          }              }
1108       }           }
1109        } /* end of parallel region */        } /* end of parallel region */
1110     } else if ( (A_block_size == 3) && (B_block_size ==3) && (C_block_size == 3) ) {     } else if ( (A_block_size == 3) && (B_block_size ==3) && (C_block_size == 3) ) {
1111        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, A_ik,B_kj, C_ij_0, C_ij_1, C_ij_2)        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, A_ik,B_kj, C_ij_0, C_ij_1, C_ij_2)
1112        {        {
1113       #pragma omp for schedule(static)           #pragma omp for schedule(static)
1114       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
1115          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
1116             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
1117             C_ij=&(C->val[ij_ptrC*C_block_size]);                 C_ij=&(C->val[ij_ptrC*C_block_size]);
1118             C_ij_0=0;                 C_ij_0=0;
1119             C_ij_1=0;                 C_ij_1=0;
1120             C_ij_2=0;                 C_ij_2=0;
1121                 ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
1122                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
1123                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
# Line 1138  void SparseMatrix_MatrixMatrixTranspose_ Line 1145  void SparseMatrix_MatrixMatrixTranspose_
1145                     kA=A->pattern->index[ik_ptrA];                     kA=A->pattern->index[ik_ptrA];
1146                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
1147                   }                   }
1148             }                 }
1149             C_ij[0]=C_ij_0;                 C_ij[0]=C_ij_0;
1150             C_ij[1]=C_ij_1;                 C_ij[1]=C_ij_1;
1151             C_ij[2]=C_ij_2;                 C_ij[2]=C_ij_2;
1152          }              }
1153       }           }
1154        } /* end of parallel region */        } /* end of parallel region */
1155     } else if ( (A_block_size == 4) && (B_block_size ==4 ) && (C_block_size == 4) ) {     } else if ( (A_block_size == 4) && (B_block_size ==4 ) && (C_block_size == 4) ) {
1156        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, A_ik,B_kj, C_ij_0, C_ij_1, C_ij_2, C_ij_3 )        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, A_ik,B_kj, C_ij_0, C_ij_1, C_ij_2, C_ij_3 )
1157        {        {
1158       #pragma omp for schedule(static)           #pragma omp for schedule(static)
1159       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
1160          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
1161             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
1162             C_ij=&(C->val[ij_ptrC*C_block_size]);                 C_ij=&(C->val[ij_ptrC*C_block_size]);
1163             C_ij_0=0;                 C_ij_0=0;
1164             C_ij_1=0;                 C_ij_1=0;
1165             C_ij_2=0;                 C_ij_2=0;
1166             C_ij_3=0;                 C_ij_3=0;
1167                 ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
1168                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
1169                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
# Line 1186  void SparseMatrix_MatrixMatrixTranspose_ Line 1193  void SparseMatrix_MatrixMatrixTranspose_
1193                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
1194                   }                   }
1195                 }                 }
1196             C_ij[0]=C_ij_0;                 C_ij[0]=C_ij_0;
1197             C_ij[1]=C_ij_1;                 C_ij[1]=C_ij_1;
1198             C_ij[2]=C_ij_2;                 C_ij[2]=C_ij_2;
1199             C_ij[3]=C_ij_3;                 C_ij[3]=C_ij_3;
1200          }              }
1201       }           }
1202        } /* end of parallel region */        } /* end of parallel region */
1203     } else {     } else {
1204        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, ib, A_ik,B_kj )        #pragma omp parallel private(C_ij, i, ij_ptrC, j, ik_ptrA, kj_ptrB, ikb, kjb, kA, kB, ib, A_ik,B_kj )
1205        {        {
1206       #pragma omp for schedule(static)           #pragma omp for schedule(static)
1207       for(i = 0; i < n; i++) {           for(i = 0; i < n; i++) {
1208          for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {              for(ij_ptrC = C->pattern->ptr[i]; ij_ptrC < C->pattern->ptr[i+1]; ++ij_ptrC) {
1209             j = C->pattern->index[ij_ptrC];                 j = C->pattern->index[ij_ptrC];
1210             C_ij=&(C->val[ij_ptrC*C_block_size]);                 C_ij=&(C->val[ij_ptrC*C_block_size]);
1211             for (ib=0; ib<C_block_size; ++ib)  C_ij[ib]=0;                 for (ib=0; ib<C_block_size; ++ib)  C_ij[ib]=0;
1212                 ik_ptrA=A->pattern->ptr[i];                 ik_ptrA=A->pattern->ptr[i];
1213                 ikb = A->pattern->ptr[i+1];                 ikb = A->pattern->ptr[i+1];
1214                 kj_ptrB=T->pattern->ptr[j];                 kj_ptrB=T->pattern->ptr[j];
# Line 1220  void SparseMatrix_MatrixMatrixTranspose_ Line 1227  void SparseMatrix_MatrixMatrixTranspose_
1227                   } else {                   } else {
1228                     A_ik=&(A->val[ik_ptrA*A_block_size]);                     A_ik=&(A->val[ik_ptrA*A_block_size]);
1229                     B_kj=&(T->val[kj_ptrB*B_block_size]);                     B_kj=&(T->val[kj_ptrB*B_block_size]);
1230             for (ib=0; ib<MIN(A_block_size, B_block_size); ++ib) C_ij[ib]+=A_ik[ib]*B_kj[ib];                     for (ib=0; ib<MIN(A_block_size, B_block_size); ++ib) C_ij[ib]+=A_ik[ib]*B_kj[ib];
1231                     ik_ptrA ++;                     ik_ptrA ++;
1232                     kj_ptrB ++;                     kj_ptrB ++;
1233                     if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;                     if (ik_ptrA >= ikb || kj_ptrB >= kjb) break;
1234                     kA=A->pattern->index[ik_ptrA];                     kA=A->pattern->index[ik_ptrA];
1235                     kB=T->pattern->index[kj_ptrB];                     kB=T->pattern->index[kj_ptrB];
1236                   }                   }
1237             }                 }
1238          }              }
1239       }           }
1240        } /* end of parallel region */        } /* end of parallel region */
1241     }     }
1242  }  }

Legend:
Removed from v.4828  
changed lines
  Added in v.4829

  ViewVC Help
Powered by ViewVC 1.1.26