/[escript]/trunk/paso/src/SparseMatrix_MatrixVector.c
ViewVC logotype

Diff of /trunk/paso/src/SparseMatrix_MatrixVector.c

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

temp/paso/src/SparseMatrix_MatrixVector.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC trunk/paso/src/SparseMatrix_MatrixVector.c revision 1565 by gross, Thu May 22 10:19:47 2008 UTC
# Line 34  void  Paso_SparseMatrix_MatrixVector_CSC Line 34  void  Paso_SparseMatrix_MatrixVector_CSC
34    
35    register index_t ir,icol,iptr,icb,irb,irow,ic;    register index_t ir,icol,iptr,icb,irb,irow,ic;
36    register double reg,reg1,reg2,reg3;    register double reg,reg1,reg2,reg3;
   #pragma omp barrier  
37    
38    if (ABS(beta)>0.) {    if (ABS(beta)>0.) {
39      if (beta != 1.) {      if (beta != 1.) {
40          #pragma omp for private(irow) schedule(static)          #pragma omp parallel for private(irow) schedule(static)
41          for (irow=0;irow < A->numRows * A->row_block_size;irow++)          for (irow=0;irow < A->numRows * A->row_block_size;irow++)
42            out[irow] *= beta;            out[irow] *= beta;
43      }      }
44    } else {    } else {
45      #pragma omp for private(irow) schedule(static)      #pragma omp parallel for private(irow) schedule(static)
46      for (irow=0;irow < A->numRows * A->row_block_size;irow++)      for (irow=0;irow < A->numRows * A->row_block_size;irow++)
47        out[irow] = 0;        out[irow] = 0;
48    }    }
# Line 52  void  Paso_SparseMatrix_MatrixVector_CSC Line 51  void  Paso_SparseMatrix_MatrixVector_CSC
51    if (ABS(alpha)>0) {    if (ABS(alpha)>0) {
52      if (A ->col_block_size==1 && A->row_block_size ==1) {      if (A ->col_block_size==1 && A->row_block_size ==1) {
53          /* TODO: parallelize (good luck!) */          /* TODO: parallelize (good luck!) */
         #pragma omp single  
54      for (icol=0;icol< A->pattern->numOutput;++icol) {      for (icol=0;icol< A->pattern->numOutput;++icol) {
55            #pragma ivdep            #pragma ivdep
56        for (iptr=A->pattern->ptr[icol];iptr<A->pattern->ptr[icol+1]; ++iptr) {        for (iptr=A->pattern->ptr[icol];iptr<A->pattern->ptr[icol+1]; ++iptr) {
# Line 61  void  Paso_SparseMatrix_MatrixVector_CSC Line 59  void  Paso_SparseMatrix_MatrixVector_CSC
59      }      }
60      } else if (A ->col_block_size==2 && A->row_block_size ==2) {      } else if (A ->col_block_size==2 && A->row_block_size ==2) {
61          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
62      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
63            #pragma ivdep            #pragma ivdep
64        for (iptr=A->pattern->ptr[ic];iptr<A->pattern->ptr[ic+1]; iptr++) {        for (iptr=A->pattern->ptr[ic];iptr<A->pattern->ptr[ic+1]; iptr++) {
65             ic=2*(A->pattern->index[iptr]);             ir=2*(A->pattern->index[iptr]);
66             out[  2*ir] += alpha * ( A->val[iptr*4  ]*in[ic] + A->val[iptr*4+2]*in[1+ic] );             out[  ir] += alpha * ( A->val[iptr*4  ]*in[ic] + A->val[iptr*4+2]*in[1+ic] );
67             out[1+2*ir] += alpha * ( A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic] );             out[1+ir] += alpha * ( A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic] );
68        }        }
69      }      }
70      } else if (A ->col_block_size==3 && A->row_block_size ==3) {      } else if (A ->col_block_size==3 && A->row_block_size ==3) {
71          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
72      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
73            #pragma ivdep            #pragma ivdep
74        for (iptr=A->pattern->ptr[ic];iptr<A->pattern->ptr[ic+1]; iptr++) {        for (iptr=A->pattern->ptr[ic];iptr<A->pattern->ptr[ic+1]; iptr++) {
75            ir=3*(A->pattern->index[iptr]);            ir=3*(A->pattern->index[iptr]);
76                out[  3*ir] += alpha * ( A->val[iptr*9  ]*in[ic] + A->val[iptr*9+3]*in[1+ic] + A->val[iptr*9+6]*in[2+ic] );                out[  ir] += alpha * ( A->val[iptr*9  ]*in[ic] + A->val[iptr*9+3]*in[1+ic] + A->val[iptr*9+6]*in[2+ic] );
77            out[1+3*ir] += alpha * ( A->val[iptr*9+1]*in[ic] + A->val[iptr*9+4]*in[1+ic] + A->val[iptr*9+7]*in[2+ic] );            out[1+ir] += alpha * ( A->val[iptr*9+1]*in[ic] + A->val[iptr*9+4]*in[1+ic] + A->val[iptr*9+7]*in[2+ic] );
78            out[2+3*ir] += alpha * ( A->val[iptr*9+2]*in[ic] + A->val[iptr*9+5]*in[1+ic] + A->val[iptr*9+8]*in[2+ic] );            out[2+ir] += alpha * ( A->val[iptr*9+2]*in[ic] + A->val[iptr*9+5]*in[1+ic] + A->val[iptr*9+8]*in[2+ic] );
79        }        }
80      }      }
81      } else {      } else {
82          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
83      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
84        for (iptr=A->pattern->ptr[ic];iptr<A->pattern->ptr[ic+1]; iptr++) {        for (iptr=A->pattern->ptr[ic];iptr<A->pattern->ptr[ic+1]; iptr++) {
85          for (irb=0;irb< A->row_block_size;irb++) {          for (irb=0;irb< A->row_block_size;irb++) {
# Line 111  void  Paso_SparseMatrix_MatrixVector_CSC Line 106  void  Paso_SparseMatrix_MatrixVector_CSC
106    
107    register index_t ir,icol,iptr,icb,irb,irow,ic;    register index_t ir,icol,iptr,icb,irb,irow,ic;
108    register double reg,reg1,reg2,reg3;    register double reg,reg1,reg2,reg3;
   #pragma omp barrier  
   
109    if (ABS(beta)>0.) {    if (ABS(beta)>0.) {
110      if (beta != 1.) {      if (beta != 1.) {
111          #pragma omp for private(irow) schedule(static)          #pragma omp parallel for private(irow) schedule(static)
112          for (irow=0;irow < A->numRows * A->row_block_size;irow++) {          for (irow=0;irow < A->numRows * A->row_block_size;irow++) {
113            out[irow] *= beta;            out[irow] *= beta;
114          }          }
115      }      }
116    } else {    } else {
117      #pragma omp for private(irow) schedule(static)      #pragma omp parallel for private(irow) schedule(static)
118      for (irow=0;irow < A->numRows * A->row_block_size;irow++)      for (irow=0;irow < A->numRows * A->row_block_size;irow++)
119        out[irow] = 0;        out[irow] = 0;
120    }    }
# Line 130  void  Paso_SparseMatrix_MatrixVector_CSC Line 123  void  Paso_SparseMatrix_MatrixVector_CSC
123    if (ABS(alpha)>0) {    if (ABS(alpha)>0) {
124      if (A ->col_block_size==1 && A->row_block_size ==1) {      if (A ->col_block_size==1 && A->row_block_size ==1) {
125          /* TODO: parallelize (good luck!) */          /* TODO: parallelize (good luck!) */
         #pragma omp single  
126      for (icol=0;icol< A->pattern->numOutput;++icol) {      for (icol=0;icol< A->pattern->numOutput;++icol) {
127            #pragma ivdep            #pragma ivdep
128        for (iptr=A->pattern->ptr[icol]-1;iptr<A->pattern->ptr[icol+1]-1; ++iptr) {        for (iptr=A->pattern->ptr[icol]-1;iptr<A->pattern->ptr[icol+1]-1; ++iptr) {
# Line 139  void  Paso_SparseMatrix_MatrixVector_CSC Line 131  void  Paso_SparseMatrix_MatrixVector_CSC
131      }      }
132      } else if (A ->col_block_size==2 && A->row_block_size ==2) {      } else if (A ->col_block_size==2 && A->row_block_size ==2) {
133          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
134      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
135        for (iptr=A->pattern->ptr[ic]-1;iptr<A->pattern->ptr[ic+1]-1; iptr++) {        for (iptr=A->pattern->ptr[ic]-1;iptr<A->pattern->ptr[ic+1]-1; iptr++) {
136             ic=2*(A->pattern->index[iptr]-1);             ir=2*(A->pattern->index[iptr]-1);
137             out[  2*ir] += alpha * ( A->val[iptr*4  ]*in[ic] + A->val[iptr*4+2]*in[1+ic] );             out[  ir] += alpha * ( A->val[iptr*4  ]*in[ic] + A->val[iptr*4+2]*in[1+ic] );
138             out[1+2*ir] += alpha * ( A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic] );             out[1+ir] += alpha * ( A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic] );
139        }        }
140      }      }
141      } else if (A ->col_block_size==3 && A->row_block_size ==3) {      } else if (A ->col_block_size==3 && A->row_block_size ==3) {
142          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
143      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
144            #pragma ivdep            #pragma ivdep
145        for (iptr=A->pattern->ptr[ic]-1;iptr<A->pattern->ptr[ic+1]-1; iptr++) {        for (iptr=A->pattern->ptr[ic]-1;iptr<A->pattern->ptr[ic+1]-1; iptr++) {
146            ir=3*(A->pattern->index[iptr]-1);            ir=3*(A->pattern->index[iptr]-1);
147                out[  3*ir] += alpha * ( A->val[iptr*9  ]*in[ic] + A->val[iptr*9+3]*in[1+ic] + A->val[iptr*9+6]*in[2+ic] );                out[  ir] += alpha * ( A->val[iptr*9  ]*in[ic] + A->val[iptr*9+3]*in[1+ic] + A->val[iptr*9+6]*in[2+ic] );
148            out[1+3*ir] += alpha * ( A->val[iptr*9+1]*in[ic] + A->val[iptr*9+4]*in[1+ic] + A->val[iptr*9+7]*in[2+ic] );            out[1+ir] += alpha * ( A->val[iptr*9+1]*in[ic] + A->val[iptr*9+4]*in[1+ic] + A->val[iptr*9+7]*in[2+ic] );
149            out[2+3*ir] += alpha * ( A->val[iptr*9+2]*in[ic] + A->val[iptr*9+5]*in[1+ic] + A->val[iptr*9+8]*in[2+ic] );            out[2+ir] += alpha * ( A->val[iptr*9+2]*in[ic] + A->val[iptr*9+5]*in[1+ic] + A->val[iptr*9+8]*in[2+ic] );
150        }        }
151      }      }
152      } else {      } else {
153          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
154      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
155        for (iptr=A->pattern->ptr[ic]-1;iptr<A->pattern->ptr[ic+1]-1; iptr++) {        for (iptr=A->pattern->ptr[ic]-1;iptr<A->pattern->ptr[ic+1]-1; iptr++) {
156          for (irb=0;irb< A->row_block_size;irb++) {          for (irb=0;irb< A->row_block_size;irb++) {
# Line 178  void  Paso_SparseMatrix_MatrixVector_CSC Line 167  void  Paso_SparseMatrix_MatrixVector_CSC
167    }    }
168    return;    return;
169  }  }
 /* CSR format with offset 0*/  
 void  Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(double alpha,  
                                                  Paso_SparseMatrix* A,  
                                                  double* in,  
                                                  double beta,  
                                                  double* out)  
 {  
     register index_t ir,icol,iptr,icb,irb,irow,ic,Aiptr;  
     register double reg,reg1,reg2,reg3,in1,in2,in3,A00,A10,A20,A01,A11,A21,A02,A12,A22;  
     if (ABS(beta)>0.) {  
       if (beta != 1.) {  
           #pragma omp for private(irow) schedule(static)  
           for (irow=0;irow < A->numRows * A->row_block_size;irow++)  
             out[irow] *= beta;  
       }  
     } else {  
       #pragma omp for private(irow) schedule(static)  
       for (irow=0;irow < A->numRows * A->row_block_size;irow++)  
         out[irow] = 0;  
     }  
     if (ABS(alpha)>0) {  
       if (A ->col_block_size==1 && A->row_block_size ==1) {  
           #pragma omp for private(irow,iptr,reg) schedule(static)  
     for (irow=0;irow< A->pattern->numOutput;++irow) {  
             reg=0.;  
             #pragma ivdep  
       for (iptr=(A->pattern->ptr[irow]);iptr<(A->pattern->ptr[irow+1]); ++iptr) {  
           reg += A->val[iptr] * in[A->pattern->index[iptr]];  
       }  
       out[irow] += alpha * reg;  
     }  
       } else if (A ->col_block_size==2 && A->row_block_size ==2) {  
           #pragma omp for private(ir,reg1,reg2,iptr,ic,Aiptr,in1,in2,A00,A10,A01,A11) schedule(static)  
     for (ir=0;ir< A->pattern->numOutput;ir++) {  
             reg1=0.;  
             reg2=0.;  
             #pragma ivdep  
       for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {  
            ic=2*(A->pattern->index[iptr]);  
                  Aiptr=iptr*4;  
                  in1=in[ic];  
                  in2=in[1+ic];  
                  A00=A->val[Aiptr  ];  
                  A10=A->val[Aiptr+1];  
                  A01=A->val[Aiptr+2];  
                  A11=A->val[Aiptr+3];  
            reg1 += A00*in1 + A01*in2;  
            reg2 += A10*in1 + A11*in2;  
       }  
       out[  2*ir] += alpha * reg1;  
       out[1+2*ir] += alpha * reg2;  
     }  
       } else if (A ->col_block_size==3 && A->row_block_size ==3) {  
           #pragma omp for private(ir,reg1,reg2,reg3,iptr,ic,Aiptr,in1,in2,in3,A00,A10,A20,A01,A11,A21,A02,A12,A22) schedule(static)  
     for (ir=0;ir< A->pattern->numOutput;ir++) {  
             reg1=0.;  
             reg2=0.;  
             reg3=0.;  
             #pragma ivdep  
       for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {  
            ic=3*(A->pattern->index[iptr]);  
                  Aiptr=iptr*9;  
                  in1=in[ic];  
                  in2=in[1+ic];  
                  in3=in[2+ic];  
                  A00=A->val[Aiptr  ];  
                  A10=A->val[Aiptr+1];  
                  A20=A->val[Aiptr+2];  
                  A01=A->val[Aiptr+3];  
                  A11=A->val[Aiptr+4];  
                  A21=A->val[Aiptr+5];  
                  A02=A->val[Aiptr+6];  
                  A12=A->val[Aiptr+7];  
                  A22=A->val[Aiptr+8];  
            reg1 += A00*in1 + A01*in2 + A02*in3;  
            reg2 += A10*in1 + A11*in2 + A12*in3;  
            reg3 += A20*in1 + A21*in2 + A22*in3;  
       }  
       out[  3*ir] += alpha * reg1;  
       out[1+3*ir] += alpha * reg2;  
       out[2+3*ir] += alpha * reg3;  
     }  
       } else {  
           #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg) schedule(static)  
     for (ir=0;ir< A->pattern->numOutput;ir++) {  
       for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {  
         for (irb=0;irb< A->row_block_size;irb++) {  
           irow=irb+A->row_block_size*ir;  
                 reg=0.;  
                 #pragma ivdep  
           for (icb=0;icb< A->col_block_size;icb++) {  
         icol=icb+A->col_block_size*(A->pattern->index[iptr]);  
         reg += A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];  
           }  
           out[irow] += alpha * reg;  
         }  
       }  
     }  
       }  
     }  
     return;  
 }  
170  /* CSR format with offset 1*/  /* CSR format with offset 1*/
171  void  Paso_SparseMatrix_MatrixVector_CSR_OFFSET1(double alpha,  void  Paso_SparseMatrix_MatrixVector_CSR_OFFSET1(double alpha,
172      Paso_SparseMatrix* A,      Paso_SparseMatrix* A,
# Line 289  void  Paso_SparseMatrix_MatrixVector_CSR Line 176  void  Paso_SparseMatrix_MatrixVector_CSR
176    
177    register index_t ir,icol,iptr,icb,irb,irow,ic;    register index_t ir,icol,iptr,icb,irb,irow,ic;
178    register double reg,reg1,reg2,reg3;    register double reg,reg1,reg2,reg3;
   #pragma omp barrier  
   
179    if (ABS(beta)>0.) {    if (ABS(beta)>0.) {
180      if (beta != 1.) {      if (beta != 1.) {
181          #pragma omp for private(irow) schedule(static)          #pragma omp parallel for private(irow) schedule(static)
182          for (irow=0;irow < A->numRows * A->row_block_size;irow++)          for (irow=0;irow < A->numRows * A->row_block_size;irow++)
183            out[irow] *= beta;            out[irow] *= beta;
184      }      }
185    } else {    } else {
186      #pragma omp for private(irow) schedule(static)      #pragma omp parallel for private(irow) schedule(static)
187      for (irow=0;irow < A->numRows * A->row_block_size;irow++)      for (irow=0;irow < A->numRows * A->row_block_size;irow++)
188        out[irow] = 0;        out[irow] = 0;
189    }    }
190    /*  do the operation: */    /*  do the operation: */
191    if (ABS(alpha)>0) {    if (ABS(alpha)>0) {
192      if (A ->col_block_size==1 && A->row_block_size ==1) {      if (A ->col_block_size==1 && A->row_block_size ==1) {
193          #pragma omp for private(irow,iptr,reg) schedule(static)          #pragma omp parallel for private(irow,iptr,reg) schedule(static)
194      for (irow=0;irow< A->pattern->numOutput;++irow) {      for (irow=0;irow< A->pattern->numOutput;++irow) {
195            reg=0.;            reg=0.;
196            #pragma ivdep            #pragma ivdep
# Line 315  void  Paso_SparseMatrix_MatrixVector_CSR Line 200  void  Paso_SparseMatrix_MatrixVector_CSR
200        out[irow] += alpha * reg;        out[irow] += alpha * reg;
201      }      }
202      } else if (A ->col_block_size==2 && A->row_block_size ==2) {      } else if (A ->col_block_size==2 && A->row_block_size ==2) {
203          #pragma omp for private(ir,reg1,reg2,iptr,ic) schedule(static)          #pragma omp parallel for private(ir,reg1,reg2,iptr,ic) schedule(static)
204      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
205            reg1=0.;            reg1=0.;
206            reg2=0.;            reg2=0.;
# Line 329  void  Paso_SparseMatrix_MatrixVector_CSR Line 214  void  Paso_SparseMatrix_MatrixVector_CSR
214        out[1+2*ir] += alpha * reg2;        out[1+2*ir] += alpha * reg2;
215      }      }
216      } else if (A ->col_block_size==3 && A->row_block_size ==3) {      } else if (A ->col_block_size==3 && A->row_block_size ==3) {
217          #pragma omp for private(ir,reg1,reg2,reg3,iptr,ic) schedule(static)          #pragma omp parallel for private(ir,reg1,reg2,reg3,iptr,ic) schedule(static)
218      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
219            reg1=0.;            reg1=0.;
220            reg2=0.;            reg2=0.;
# Line 346  void  Paso_SparseMatrix_MatrixVector_CSR Line 231  void  Paso_SparseMatrix_MatrixVector_CSR
231        out[2+3*ir] += alpha * reg3;        out[2+3*ir] += alpha * reg3;
232      }      }
233      } else {      } else {
234          #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg) schedule(static)          #pragma omp parallel for private(ir,iptr,irb,icb,irow,icol,reg) schedule(static)
235      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
236        for (iptr=A->pattern->ptr[ir]-1;iptr<A->pattern->ptr[ir+1]-1; iptr++) {        for (iptr=A->pattern->ptr[ir]-1;iptr<A->pattern->ptr[ir+1]-1; iptr++) {
237          for (irb=0;irb< A->row_block_size;irb++) {          for (irb=0;irb< A->row_block_size;irb++) {
# Line 365  void  Paso_SparseMatrix_MatrixVector_CSR Line 250  void  Paso_SparseMatrix_MatrixVector_CSR
250    }    }
251    return;    return;
252  }  }
253    /* CSR format with offset 0*/
254    void  Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(double alpha,
255                                                     Paso_SparseMatrix* A,
256                                                     double* in,
257                                                     double beta,
258                                                     double* out)
259    {
260    #define PASO_DYNAMIC_SCHEDULING_MVM
261    
262        char* chksz_chr=NULL;
263        dim_t chunk_size=-1;
264        dim_t nrow=A->numRows;
265        dim_t np, len, rest, irow, local_n, p, n_chunks;
266        np=omp_get_max_threads();
267    #if defined PASO_DYNAMIC_SCHEDULING_MVM && defined __OPENMP
268        chksz_chr=getenv("PASO_CHUNK_SIZE_MVM");
269        if (chksz_chr!=NULL) sscanf(chksz_chr, "%d",&chunk_size);
270    #endif
271        
272       if (chunk_size<1 || np <=1) {
273            #pragma omp parallel private(irow, len, rest, local_n)
274            {
275               len=nrow/np;
276               rest=nrow-len*np;
277               #pragma omp for private(p) schedule(static)
278               for (p=0; p<np;p++) {
279                  irow=len*p+MIN(p,rest);
280                  local_n=len+(p<rest ? 1 :0 );
281                  Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_stripe(alpha,
282                                                                    local_n,
283                                                                    A->row_block_size,
284                                                                    A->col_block_size,
285                                                                    &(A->pattern->ptr[irow]),
286                                                                    A->pattern->index, A->val, in, beta, &out[irow*A->row_block_size]);
287               }
288            }
289       } else {
290          #pragma omp parallel private(n_chunks,irow,local_n)
291          {
292             n_chunks=nrow/chunk_size;
293             if (n_chunks*chunk_size<nrow) n_chunks+=1;
294             #pragma omp for private(p) schedule(dynamic,1)
295             for (p=0; p<n_chunks;p++) {
296               irow=chunk_size*p;
297               local_n=MIN(chunk_size,nrow-chunk_size*p);
298               Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_stripe(alpha,
299                                                                 local_n,
300                                                                 A->row_block_size,
301                                                                 A->col_block_size,
302                                                                 &(A->pattern->ptr[irow]),
303                                                                 A->pattern->index, A->val, in, beta, &out[irow*A->row_block_size]);
304            }
305          }
306       }
307    }
308    /* CSR format with offset 0*/
309    void  Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_stripe(double alpha,
310                                                            dim_t nRows,
311                                                            dim_t row_block_size,
312                                                            dim_t col_block_size,
313                                                            index_t* ptr,
314                                                            index_t* index,
315                                                            double* val,
316                                                            double* in,
317                                                            double beta,
318                                                            double* out)
319    {
320        register index_t ir,icol,iptr,icb,irb,irow,ic,Aiptr;
321        register double reg,reg1,reg2,reg3,in1,in2,in3,A00,A10,A20,A01,A11,A21,A02,A12,A22;
322        dim_t block_size;
323        if (ABS(beta)>0.) {
324          if (beta != 1.) {
325              for (irow=0;irow < nRows * row_block_size;irow++)
326                out[irow] *= beta;
327          }
328        } else {
329          for (irow=0;irow < nRows * row_block_size;irow++)
330            out[irow] = 0;
331        }
332        if (ABS(alpha)>0) {
333          if (col_block_size==1 && row_block_size ==1) {
334        for (irow=0;irow< nRows;++irow) {
335              reg=0.;
336              #pragma ivdep
337          for (iptr=ptr[irow];iptr<ptr[irow+1]; ++iptr) {
338              reg += val[iptr] * in[index[iptr]];
339          }
340          out[irow] += alpha * reg;
341        }
342          } else if (col_block_size==2 && row_block_size ==2) {
343        for (ir=0;ir< nRows;ir++) {
344              reg1=0.;
345              reg2=0.;
346              #pragma ivdep
347          for (iptr=ptr[ir];iptr<ptr[ir+1]; iptr++) {
348               ic=2*index[iptr];
349                     Aiptr=iptr*4;
350                     in1=in[ic];
351                     in2=in[1+ic];
352                     A00=val[Aiptr  ];
353                     A10=val[Aiptr+1];
354                     A01=val[Aiptr+2];
355                     A11=val[Aiptr+3];
356               reg1 += A00*in1 + A01*in2;
357               reg2 += A10*in1 + A11*in2;
358          }
359          out[  2*ir] += alpha * reg1;
360          out[1+2*ir] += alpha * reg2;
361        }
362          } else if (col_block_size==3 && row_block_size ==3) {
363        for (ir=0;ir< nRows;ir++) {
364              reg1=0.;
365              reg2=0.;
366              reg3=0.;
367              #pragma ivdep
368          for (iptr=ptr[ir];iptr<ptr[ir+1]; iptr++) {
369               ic=3*index[iptr];
370                     Aiptr=iptr*9;
371                     in1=in[ic];
372                     in2=in[1+ic];
373                     in3=in[2+ic];
374                     A00=val[Aiptr  ];
375                     A10=val[Aiptr+1];
376                     A20=val[Aiptr+2];
377                     A01=val[Aiptr+3];
378                     A11=val[Aiptr+4];
379                     A21=val[Aiptr+5];
380                     A02=val[Aiptr+6];
381                     A12=val[Aiptr+7];
382                     A22=val[Aiptr+8];
383               reg1 += A00*in1 + A01*in2 + A02*in3;
384               reg2 += A10*in1 + A11*in2 + A12*in3;
385               reg3 += A20*in1 + A21*in2 + A22*in3;
386          }
387          out[  3*ir] += alpha * reg1;
388          out[1+3*ir] += alpha * reg2;
389          out[2+3*ir] += alpha * reg3;
390        }
391          } else {
392            block_size=col_block_size*row_block_size;
393        for (ir=0;ir< nRows;ir++) {
394          for (iptr=ptr[ir];iptr<ptr[ir+1]; iptr++) {
395            for (irb=0;irb< row_block_size;irb++) {
396              irow=irb+row_block_size*ir;
397                  reg=0.;
398                  #pragma ivdep
399              for (icb=0;icb< col_block_size;icb++) {
400            icol=icb+col_block_size*index[iptr];
401            reg += val[iptr*block_size+irb+row_block_size*icb] * in[icol];
402              }
403              out[irow] += alpha * reg;
404            }
405          }
406        }
407          }
408        }
409        return;
410    }

Legend:
Removed from v.1387  
changed lines
  Added in v.1565

  ViewVC Help
Powered by ViewVC 1.1.26