/[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 1639 by gross, Mon Jul 14 08:55:25 2008 UTC
# Line 24  Line 24 
24  /**************************************************************/  /**************************************************************/
25    
26  #include "SparseMatrix.h"  #include "SparseMatrix.h"
27    #ifdef _OPENMP
28    #include <omp.h>
29    #endif
30    
31  /* CSC format with offset 0*/  /* CSC format with offset 0*/
32  void  Paso_SparseMatrix_MatrixVector_CSC_OFFSET0(double alpha,  void  Paso_SparseMatrix_MatrixVector_CSC_OFFSET0(const double alpha,
33                                                   Paso_SparseMatrix* A,                                                   const Paso_SparseMatrix* A,
34                                                   double* in,                                                   const double* in,
35                                                   double beta,                                                   const double beta,
36                                                   double* out) {                                                   double* out) {
37    
38    register index_t ir,icol,iptr,icb,irb,irow,ic;    register index_t ir,icol,iptr,icb,irb,irow,ic;
   register double reg,reg1,reg2,reg3;  
   #pragma omp barrier  
39    
40    if (ABS(beta)>0.) {    if (ABS(beta)>0.) {
41      if (beta != 1.) {      if (beta != 1.) {
42          #pragma omp for private(irow) schedule(static)          #pragma omp parallel for private(irow) schedule(static)
43          for (irow=0;irow < A->numRows * A->row_block_size;irow++)          for (irow=0;irow < A->numRows * A->row_block_size;irow++)
44            out[irow] *= beta;            out[irow] *= beta;
45      }      }
46    } else {    } else {
47      #pragma omp for private(irow) schedule(static)      #pragma omp parallel for private(irow) schedule(static)
48      for (irow=0;irow < A->numRows * A->row_block_size;irow++)      for (irow=0;irow < A->numRows * A->row_block_size;irow++)
49        out[irow] = 0;        out[irow] = 0;
50    }    }
# Line 52  void  Paso_SparseMatrix_MatrixVector_CSC Line 53  void  Paso_SparseMatrix_MatrixVector_CSC
53    if (ABS(alpha)>0) {    if (ABS(alpha)>0) {
54      if (A ->col_block_size==1 && A->row_block_size ==1) {      if (A ->col_block_size==1 && A->row_block_size ==1) {
55          /* TODO: parallelize (good luck!) */          /* TODO: parallelize (good luck!) */
         #pragma omp single  
56      for (icol=0;icol< A->pattern->numOutput;++icol) {      for (icol=0;icol< A->pattern->numOutput;++icol) {
57            #pragma ivdep            #pragma ivdep
58        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 61  void  Paso_SparseMatrix_MatrixVector_CSC
61      }      }
62      } else if (A ->col_block_size==2 && A->row_block_size ==2) {      } else if (A ->col_block_size==2 && A->row_block_size ==2) {
63          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
64      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
65            #pragma ivdep            #pragma ivdep
66        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++) {
67             ic=2*(A->pattern->index[iptr]);             ir=2*(A->pattern->index[iptr]);
68             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] );
69             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] );
70        }        }
71      }      }
72      } else if (A ->col_block_size==3 && A->row_block_size ==3) {      } else if (A ->col_block_size==3 && A->row_block_size ==3) {
73          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
74      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
75            #pragma ivdep            #pragma ivdep
76        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++) {
77            ir=3*(A->pattern->index[iptr]);            ir=3*(A->pattern->index[iptr]);
78                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] );
79            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] );
80            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] );
81        }        }
82      }      }
83      } else {      } else {
84          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
85      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
86        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++) {
87          for (irb=0;irb< A->row_block_size;irb++) {          for (irb=0;irb< A->row_block_size;irb++) {
# Line 103  void  Paso_SparseMatrix_MatrixVector_CSC Line 100  void  Paso_SparseMatrix_MatrixVector_CSC
100  }  }
101    
102  /* CSC format with offset 1*/  /* CSC format with offset 1*/
103  void  Paso_SparseMatrix_MatrixVector_CSC_OFFSET1(double alpha,  void  Paso_SparseMatrix_MatrixVector_CSC_OFFSET1(const double alpha,
104                                                   Paso_SparseMatrix* A,                                                   const Paso_SparseMatrix* A,
105                                                   double* in,                                                   const double* in,
106                                                   double beta,                                                   const double beta,
107                                                   double* out) {                                                   double* out) {
108    
109    register index_t ir,icol,iptr,icb,irb,irow,ic;    register index_t ir,icol,iptr,icb,irb,irow,ic;
   register double reg,reg1,reg2,reg3;  
   #pragma omp barrier  
110    
111    if (ABS(beta)>0.) {    if (ABS(beta)>0.) {
112      if (beta != 1.) {      if (beta != 1.) {
113          #pragma omp for private(irow) schedule(static)          #pragma omp parallel for private(irow) schedule(static)
114          for (irow=0;irow < A->numRows * A->row_block_size;irow++) {          for (irow=0;irow < A->numRows * A->row_block_size;irow++) {
115            out[irow] *= beta;            out[irow] *= beta;
116          }          }
117      }      }
118    } else {    } else {
119      #pragma omp for private(irow) schedule(static)      #pragma omp parallel for private(irow) schedule(static)
120      for (irow=0;irow < A->numRows * A->row_block_size;irow++)      for (irow=0;irow < A->numRows * A->row_block_size;irow++)
121        out[irow] = 0;        out[irow] = 0;
122    }    }
# Line 130  void  Paso_SparseMatrix_MatrixVector_CSC Line 125  void  Paso_SparseMatrix_MatrixVector_CSC
125    if (ABS(alpha)>0) {    if (ABS(alpha)>0) {
126      if (A ->col_block_size==1 && A->row_block_size ==1) {      if (A ->col_block_size==1 && A->row_block_size ==1) {
127          /* TODO: parallelize (good luck!) */          /* TODO: parallelize (good luck!) */
         #pragma omp single  
128      for (icol=0;icol< A->pattern->numOutput;++icol) {      for (icol=0;icol< A->pattern->numOutput;++icol) {
129            #pragma ivdep            #pragma ivdep
130        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 133  void  Paso_SparseMatrix_MatrixVector_CSC
133      }      }
134      } else if (A ->col_block_size==2 && A->row_block_size ==2) {      } else if (A ->col_block_size==2 && A->row_block_size ==2) {
135          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
136      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
137        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++) {
138             ic=2*(A->pattern->index[iptr]-1);             ir=2*(A->pattern->index[iptr]-1);
139             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] );
140             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] );
141        }        }
142      }      }
143      } else if (A ->col_block_size==3 && A->row_block_size ==3) {      } else if (A ->col_block_size==3 && A->row_block_size ==3) {
144          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
145      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
146            #pragma ivdep            #pragma ivdep
147        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++) {
148            ir=3*(A->pattern->index[iptr]-1);            ir=3*(A->pattern->index[iptr]-1);
149                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] );
150            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] );
151            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] );
152        }        }
153      }      }
154      } else {      } else {
155          /* TODO: parallelize */          /* TODO: parallelize */
         #pragma omp single  
156      for (ic=0;ic< A->pattern->numOutput;ic++) {      for (ic=0;ic< A->pattern->numOutput;ic++) {
157        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++) {
158          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 169  void  Paso_SparseMatrix_MatrixVector_CSC
169    }    }
170    return;    return;
171  }  }
 /* 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;  
 }  
172  /* CSR format with offset 1*/  /* CSR format with offset 1*/
173  void  Paso_SparseMatrix_MatrixVector_CSR_OFFSET1(double alpha,  void  Paso_SparseMatrix_MatrixVector_CSR_OFFSET1(const double alpha,
174      Paso_SparseMatrix* A,      const Paso_SparseMatrix* A,
175      double* in,      const double* in,
176      double beta,      const double beta,
177      double* out) {      double* out) {
178    
179    register index_t ir,icol,iptr,icb,irb,irow,ic;    register index_t ir,icol,iptr,icb,irb,irow,ic;
180    register double reg,reg1,reg2,reg3;    register double reg,reg1,reg2,reg3;
   #pragma omp barrier  
   
181    if (ABS(beta)>0.) {    if (ABS(beta)>0.) {
182      if (beta != 1.) {      if (beta != 1.) {
183          #pragma omp for private(irow) schedule(static)          #pragma omp parallel for private(irow) schedule(static)
184          for (irow=0;irow < A->numRows * A->row_block_size;irow++)          for (irow=0;irow < A->numRows * A->row_block_size;irow++)
185            out[irow] *= beta;            out[irow] *= beta;
186      }      }
187    } else {    } else {
188      #pragma omp for private(irow) schedule(static)      #pragma omp parallel for private(irow) schedule(static)
189      for (irow=0;irow < A->numRows * A->row_block_size;irow++)      for (irow=0;irow < A->numRows * A->row_block_size;irow++)
190        out[irow] = 0;        out[irow] = 0;
191    }    }
192    /*  do the operation: */    /*  do the operation: */
193    if (ABS(alpha)>0) {    if (ABS(alpha)>0) {
194      if (A ->col_block_size==1 && A->row_block_size ==1) {      if (A ->col_block_size==1 && A->row_block_size ==1) {
195          #pragma omp for private(irow,iptr,reg) schedule(static)          #pragma omp parallel for private(irow,iptr,reg) schedule(static)
196      for (irow=0;irow< A->pattern->numOutput;++irow) {      for (irow=0;irow< A->pattern->numOutput;++irow) {
197            reg=0.;            reg=0.;
198            #pragma ivdep            #pragma ivdep
# Line 315  void  Paso_SparseMatrix_MatrixVector_CSR Line 202  void  Paso_SparseMatrix_MatrixVector_CSR
202        out[irow] += alpha * reg;        out[irow] += alpha * reg;
203      }      }
204      } else if (A ->col_block_size==2 && A->row_block_size ==2) {      } else if (A ->col_block_size==2 && A->row_block_size ==2) {
205          #pragma omp for private(ir,reg1,reg2,iptr,ic) schedule(static)          #pragma omp parallel for private(ir,reg1,reg2,iptr,ic) schedule(static)
206      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
207            reg1=0.;            reg1=0.;
208            reg2=0.;            reg2=0.;
# Line 329  void  Paso_SparseMatrix_MatrixVector_CSR Line 216  void  Paso_SparseMatrix_MatrixVector_CSR
216        out[1+2*ir] += alpha * reg2;        out[1+2*ir] += alpha * reg2;
217      }      }
218      } else if (A ->col_block_size==3 && A->row_block_size ==3) {      } else if (A ->col_block_size==3 && A->row_block_size ==3) {
219          #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)
220      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
221            reg1=0.;            reg1=0.;
222            reg2=0.;            reg2=0.;
# Line 346  void  Paso_SparseMatrix_MatrixVector_CSR Line 233  void  Paso_SparseMatrix_MatrixVector_CSR
233        out[2+3*ir] += alpha * reg3;        out[2+3*ir] += alpha * reg3;
234      }      }
235      } else {      } else {
236          #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)
237      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
238        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++) {
239          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 252  void  Paso_SparseMatrix_MatrixVector_CSR
252    }    }
253    return;    return;
254  }  }
255    /* CSR format with offset 0*/
256    void  Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(const double alpha,
257                                                     const Paso_SparseMatrix* A,
258                                                     const double* in,
259                                                     const double beta,
260                                                     double* out)
261    {
262    /*#define PASO_DYNAMIC_SCHEDULING_MVM */
263    
264    #if defined PASO_DYNAMIC_SCHEDULING_MVM && defined _OPENMP
265    #define USE_DYNAMIC_SCHEDULING
266    #endif
267    
268        char* chksz_chr=NULL;
269        dim_t chunk_size=1;
270        dim_t nrow=A->numRows;
271        dim_t np=1, len, rest, irow, local_n, p, n_chunks;
272    #ifdef _OPENMP
273        np=omp_get_max_threads();
274    #endif
275    #ifdef USE_DYNAMIC_SCHEDULING
276        chksz_chr=getenv("PASO_CHUNK_SIZE_MVM");
277        if (chksz_chr!=NULL) sscanf(chksz_chr, "%d",&chunk_size);
278        chunk_size=MIN(MAX(1,chunk_size),nrow/np);
279        n_chunks=nrow/chunk_size;
280        if (n_chunks*chunk_size<nrow) n_chunks+=1;
281    #else
282        len=nrow/np;
283        rest=nrow-len*np;
284    #endif
285        
286         #pragma omp parallel private(irow, p, local_n)
287         {
288            #ifdef USE_DYNAMIC_SCHEDULING
289              #pragma omp for schedule(dynamic,1)
290              for (p=0; p<n_chunks;p++) {
291                irow=chunk_size*p;
292                local_n=MIN(chunk_size,nrow-chunk_size*p);
293            #else
294                #pragma omp for schedule(static)
295                for (p=0; p<np;p++) {
296                   irow=len*p+MIN(p,rest);
297                   local_n=len+(p<rest ? 1 :0 );
298            #endif
299            Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_stripe(alpha,
300                                                              local_n,
301                                                              A->row_block_size,
302                                                              A->col_block_size,
303                                                              &(A->pattern->ptr[irow]),
304                                                               A->pattern->index, A->val, in, beta, &out[irow*A->row_block_size]);
305            #ifdef USE_DYNAMIC_SCHEDULING
306              }
307            #else
308              }
309            #endif
310       }
311    }
312    /* CSR format with offset 0*/
313    void  Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_stripe(const double alpha,
314                                                            const dim_t nRows,
315                                                            const dim_t row_block_size,
316                                                            const dim_t col_block_size,
317                                                            const index_t* ptr,
318                                                            const index_t* index,
319                                                            const double* val,
320                                                            const double* in,
321                                                            const double beta,
322                                                            double* out)
323    {
324        register index_t ir,icol,iptr,icb,irb,irow,ic,Aiptr;
325        register double reg,reg1,reg2,reg3,in1,in2,in3,A00,A10,A20,A01,A11,A21,A02,A12,A22;
326        dim_t block_size;
327        if (ABS(beta)>0.) {
328          if (beta != 1.) {
329              for (irow=0;irow < nRows * row_block_size;irow++)
330                out[irow] *= beta;
331          }
332        } else {
333          for (irow=0;irow < nRows * row_block_size;irow++)
334            out[irow] = 0;
335        }
336        if (ABS(alpha)>0) {
337          if (col_block_size==1 && row_block_size ==1) {
338        for (irow=0;irow< nRows;++irow) {
339              reg=0.;
340              #pragma ivdep
341          for (iptr=ptr[irow];iptr<ptr[irow+1]; ++iptr) {
342              reg += val[iptr] * in[index[iptr]];
343          }
344          out[irow] += alpha * reg;
345        }
346          } else if (col_block_size==2 && row_block_size ==2) {
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        for (ir=0;ir< nRows;ir++) {
368              reg1=0.;
369              reg2=0.;
370              reg3=0.;
371              #pragma ivdep
372          for (iptr=ptr[ir];iptr<ptr[ir+1]; iptr++) {
373               ic=3*index[iptr];
374                     Aiptr=iptr*9;
375                     in1=in[ic];
376                     in2=in[1+ic];
377                     in3=in[2+ic];
378                     A00=val[Aiptr  ];
379                     A10=val[Aiptr+1];
380                     A20=val[Aiptr+2];
381                     A01=val[Aiptr+3];
382                     A11=val[Aiptr+4];
383                     A21=val[Aiptr+5];
384                     A02=val[Aiptr+6];
385                     A12=val[Aiptr+7];
386                     A22=val[Aiptr+8];
387               reg1 += A00*in1 + A01*in2 + A02*in3;
388               reg2 += A10*in1 + A11*in2 + A12*in3;
389               reg3 += A20*in1 + A21*in2 + A22*in3;
390          }
391          out[  3*ir] += alpha * reg1;
392          out[1+3*ir] += alpha * reg2;
393          out[2+3*ir] += alpha * reg3;
394        }
395          } else {
396            block_size=col_block_size*row_block_size;
397        for (ir=0;ir< nRows;ir++) {
398          for (iptr=ptr[ir];iptr<ptr[ir+1]; iptr++) {
399            for (irb=0;irb< row_block_size;irb++) {
400              irow=irb+row_block_size*ir;
401                  reg=0.;
402                  #pragma ivdep
403              for (icb=0;icb< col_block_size;icb++) {
404            icol=icb+col_block_size*index[iptr];
405            reg += val[iptr*block_size+irb+row_block_size*icb] * in[icol];
406              }
407              out[irow] += alpha * reg;
408            }
409          }
410        }
411          }
412        }
413        return;
414    }

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

  ViewVC Help
Powered by ViewVC 1.1.26