/[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 1511 by gross, Mon Apr 14 23:09:38 2008 UTC revision 1556 by gross, Mon May 12 00:54:58 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++) {
# Line 72  void  Paso_SparseMatrix_MatrixVector_CSC Line 69  void  Paso_SparseMatrix_MatrixVector_CSC
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++) {
# Line 84  void  Paso_SparseMatrix_MatrixVector_CSC Line 80  void  Paso_SparseMatrix_MatrixVector_CSC
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             ir=2*(A->pattern->index[iptr]-1);             ir=2*(A->pattern->index[iptr]-1);
# Line 149  void  Paso_SparseMatrix_MatrixVector_CSC Line 140  void  Paso_SparseMatrix_MatrixVector_CSC
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++) {
# Line 161  void  Paso_SparseMatrix_MatrixVector_CSC Line 151  void  Paso_SparseMatrix_MatrixVector_CSC
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 189  void  Paso_SparseMatrix_MatrixVector_CSR Line 178  void  Paso_SparseMatrix_MatrixVector_CSR
178      register double reg,reg1,reg2,reg3,in1,in2,in3,A00,A10,A20,A01,A11,A21,A02,A12,A22;      register double reg,reg1,reg2,reg3,in1,in2,in3,A00,A10,A20,A01,A11,A21,A02,A12,A22;
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      if (ABS(alpha)>0) {      if (ABS(alpha)>0) {
191        if (A ->col_block_size==1 && A->row_block_size ==1) {        if (A ->col_block_size==1 && A->row_block_size ==1) {
192            #pragma omp for private(irow,iptr,reg) schedule(static)            #pragma omp parallel for private(irow,iptr,reg) schedule(static)
193      for (irow=0;irow< A->pattern->numOutput;++irow) {      for (irow=0;irow< A->pattern->numOutput;++irow) {
194              reg=0.;              reg=0.;
195              #pragma ivdep              #pragma ivdep
# Line 210  void  Paso_SparseMatrix_MatrixVector_CSR Line 199  void  Paso_SparseMatrix_MatrixVector_CSR
199        out[irow] += alpha * reg;        out[irow] += alpha * reg;
200      }      }
201        } else if (A ->col_block_size==2 && A->row_block_size ==2) {        } else if (A ->col_block_size==2 && A->row_block_size ==2) {
202            #pragma omp for private(ir,reg1,reg2,iptr,ic,Aiptr,in1,in2,A00,A10,A01,A11) schedule(static)          #pragma omp parallel for private(ir,reg1,reg2,iptr,ic,Aiptr,in1,in2,A00,A10,A01,A11) schedule(static)
203      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
204              reg1=0.;            reg1=0.;
205              reg2=0.;            reg2=0.;
206              #pragma ivdep            #pragma ivdep
207        for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {        for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {
208             ic=2*(A->pattern->index[iptr]);             ic=2*(A->pattern->index[iptr]);
209                   Aiptr=iptr*4;                   Aiptr=iptr*4;
# Line 231  void  Paso_SparseMatrix_MatrixVector_CSR Line 220  void  Paso_SparseMatrix_MatrixVector_CSR
220        out[1+2*ir] += alpha * reg2;        out[1+2*ir] += alpha * reg2;
221      }      }
222        } else if (A ->col_block_size==3 && A->row_block_size ==3) {        } else if (A ->col_block_size==3 && A->row_block_size ==3) {
223            #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)          #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)
224      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
225              reg1=0.;            reg1=0.;
226              reg2=0.;            reg2=0.;
227              reg3=0.;            reg3=0.;
228              #pragma ivdep            #pragma ivdep
229        for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {        for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {
230             ic=3*(A->pattern->index[iptr]);             ic=3*(A->pattern->index[iptr]);
231                   Aiptr=iptr*9;                   Aiptr=iptr*9;
# Line 261  void  Paso_SparseMatrix_MatrixVector_CSR Line 250  void  Paso_SparseMatrix_MatrixVector_CSR
250        out[2+3*ir] += alpha * reg3;        out[2+3*ir] += alpha * reg3;
251      }      }
252        } else {        } else {
253            #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)
254      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
255        for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {        for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {
256          for (irb=0;irb< A->row_block_size;irb++) {          for (irb=0;irb< A->row_block_size;irb++) {
257            irow=irb+A->row_block_size*ir;            irow=irb+A->row_block_size*ir;
258                  reg=0.;                reg=0.;
259                  #pragma ivdep                #pragma ivdep
260            for (icb=0;icb< A->col_block_size;icb++) {            for (icb=0;icb< A->col_block_size;icb++) {
261          icol=icb+A->col_block_size*(A->pattern->index[iptr]);          icol=icb+A->col_block_size*(A->pattern->index[iptr]);
262          reg += A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];          reg += A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];
# Line 289  void  Paso_SparseMatrix_MatrixVector_CSR Line 278  void  Paso_SparseMatrix_MatrixVector_CSR
278    
279    register index_t ir,icol,iptr,icb,irb,irow,ic;    register index_t ir,icol,iptr,icb,irb,irow,ic;
280    register double reg,reg1,reg2,reg3;    register double reg,reg1,reg2,reg3;
   #pragma omp barrier  
   
281    if (ABS(beta)>0.) {    if (ABS(beta)>0.) {
282      if (beta != 1.) {      if (beta != 1.) {
283          #pragma omp for private(irow) schedule(static)          #pragma omp parallel for private(irow) schedule(static)
284          for (irow=0;irow < A->numRows * A->row_block_size;irow++)          for (irow=0;irow < A->numRows * A->row_block_size;irow++)
285            out[irow] *= beta;            out[irow] *= beta;
286      }      }
287    } else {    } else {
288      #pragma omp for private(irow) schedule(static)      #pragma omp parallel for private(irow) schedule(static)
289      for (irow=0;irow < A->numRows * A->row_block_size;irow++)      for (irow=0;irow < A->numRows * A->row_block_size;irow++)
290        out[irow] = 0;        out[irow] = 0;
291    }    }
292    /*  do the operation: */    /*  do the operation: */
293    if (ABS(alpha)>0) {    if (ABS(alpha)>0) {
294      if (A ->col_block_size==1 && A->row_block_size ==1) {      if (A ->col_block_size==1 && A->row_block_size ==1) {
295          #pragma omp for private(irow,iptr,reg) schedule(static)          #pragma omp parallel for private(irow,iptr,reg) schedule(static)
296      for (irow=0;irow< A->pattern->numOutput;++irow) {      for (irow=0;irow< A->pattern->numOutput;++irow) {
297            reg=0.;            reg=0.;
298            #pragma ivdep            #pragma ivdep
# Line 315  void  Paso_SparseMatrix_MatrixVector_CSR Line 302  void  Paso_SparseMatrix_MatrixVector_CSR
302        out[irow] += alpha * reg;        out[irow] += alpha * reg;
303      }      }
304      } else if (A ->col_block_size==2 && A->row_block_size ==2) {      } else if (A ->col_block_size==2 && A->row_block_size ==2) {
305          #pragma omp for private(ir,reg1,reg2,iptr,ic) schedule(static)          #pragma omp parallel for private(ir,reg1,reg2,iptr,ic) schedule(static)
306      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
307            reg1=0.;            reg1=0.;
308            reg2=0.;            reg2=0.;
# Line 329  void  Paso_SparseMatrix_MatrixVector_CSR Line 316  void  Paso_SparseMatrix_MatrixVector_CSR
316        out[1+2*ir] += alpha * reg2;        out[1+2*ir] += alpha * reg2;
317      }      }
318      } else if (A ->col_block_size==3 && A->row_block_size ==3) {      } else if (A ->col_block_size==3 && A->row_block_size ==3) {
319          #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)
320      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
321            reg1=0.;            reg1=0.;
322            reg2=0.;            reg2=0.;
# Line 346  void  Paso_SparseMatrix_MatrixVector_CSR Line 333  void  Paso_SparseMatrix_MatrixVector_CSR
333        out[2+3*ir] += alpha * reg3;        out[2+3*ir] += alpha * reg3;
334      }      }
335      } else {      } else {
336          #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)
337      for (ir=0;ir< A->pattern->numOutput;ir++) {      for (ir=0;ir< A->pattern->numOutput;ir++) {
338        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++) {
339          for (irb=0;irb< A->row_block_size;irb++) {          for (irb=0;irb< A->row_block_size;irb++) {

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

  ViewVC Help
Powered by ViewVC 1.1.26