/[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

revision 1556 by gross, Mon May 12 00:54:58 2008 UTC revision 1564 by gross, Thu May 22 09:31:33 2008 UTC
# Line 167  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 parallel for private(irow) schedule(static)  
           for (irow=0;irow < A->numRows * A->row_block_size;irow++)  
             out[irow] *= beta;  
       }  
     } else {  
       #pragma omp parallel 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 parallel 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 parallel 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 parallel 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 parallel 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 352  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              #pragma omp for private(irow) schedule(static)
326              for (irow=0;irow < nRows * row_block_size;irow++)
327                out[irow] *= beta;
328          }
329        } else {
330          #pragma omp for private(irow) schedule(static)
331          for (irow=0;irow < nRows * row_block_size;irow++)
332            out[irow] = 0;
333        }
334        if (ABS(alpha)>0) {
335          if (col_block_size==1 && row_block_size ==1) {
336              #pragma omp for private(irow,iptr,reg) schedule(static)
337        for (irow=0;irow< nRows;++irow) {
338              reg=0.;
339              #pragma ivdep
340          for (iptr=ptr[irow];iptr<ptr[irow+1]; ++iptr) {
341              reg += val[iptr] * in[index[iptr]];
342          }
343          out[irow] += alpha * reg;
344        }
345          } else if (col_block_size==2 && row_block_size ==2) {
346            #pragma omp for private(ir,reg1,reg2,iptr,ic,Aiptr,in1,in2,A00,A10,A01,A11) schedule(static)
347        for (ir=0;ir< nRows;ir++) {
348              reg1=0.;
349              reg2=0.;
350              #pragma ivdep
351          for (iptr=ptr[ir];iptr<ptr[ir+1]; iptr++) {
352               ic=2*index[iptr];
353                     Aiptr=iptr*4;
354                     in1=in[ic];
355                     in2=in[1+ic];
356                     A00=val[Aiptr  ];
357                     A10=val[Aiptr+1];
358                     A01=val[Aiptr+2];
359                     A11=val[Aiptr+3];
360               reg1 += A00*in1 + A01*in2;
361               reg2 += A10*in1 + A11*in2;
362          }
363          out[  2*ir] += alpha * reg1;
364          out[1+2*ir] += alpha * reg2;
365        }
366          } else if (col_block_size==3 && row_block_size ==3) {
367            #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)
368        for (ir=0;ir< nRows;ir++) {
369              reg1=0.;
370              reg2=0.;
371              reg3=0.;
372              #pragma ivdep
373          for (iptr=ptr[ir];iptr<ptr[ir+1]; iptr++) {
374               ic=3*index[iptr];
375                     Aiptr=iptr*9;
376                     in1=in[ic];
377                     in2=in[1+ic];
378                     in3=in[2+ic];
379                     A00=val[Aiptr  ];
380                     A10=val[Aiptr+1];
381                     A20=val[Aiptr+2];
382                     A01=val[Aiptr+3];
383                     A11=val[Aiptr+4];
384                     A21=val[Aiptr+5];
385                     A02=val[Aiptr+6];
386                     A12=val[Aiptr+7];
387                     A22=val[Aiptr+8];
388               reg1 += A00*in1 + A01*in2 + A02*in3;
389               reg2 += A10*in1 + A11*in2 + A12*in3;
390               reg3 += A20*in1 + A21*in2 + A22*in3;
391          }
392          out[  3*ir] += alpha * reg1;
393          out[1+3*ir] += alpha * reg2;
394          out[2+3*ir] += alpha * reg3;
395        }
396          } else {
397            block_size=col_block_size*row_block_size;
398            #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg) schedule(static)
399        for (ir=0;ir< nRows;ir++) {
400          for (iptr=ptr[ir];iptr<ptr[ir+1]; iptr++) {
401            for (irb=0;irb< row_block_size;irb++) {
402              irow=irb+row_block_size*ir;
403                  reg=0.;
404                  #pragma ivdep
405              for (icb=0;icb< col_block_size;icb++) {
406            icol=icb+col_block_size*index[iptr];
407            reg += val[iptr*block_size+irb+row_block_size*icb] * in[icol];
408              }
409              out[irow] += alpha * reg;
410            }
411          }
412        }
413          }
414        }
415        return;
416    }

Legend:
Removed from v.1556  
changed lines
  Added in v.1564

  ViewVC Help
Powered by ViewVC 1.1.26