/[escript]/trunk/esys2/finley/src/finleyC/System_MVM.c
ViewVC logotype

Diff of /trunk/esys2/finley/src/finleyC/System_MVM.c

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

revision 112 by jgs, Wed Dec 15 07:08:39 2004 UTC revision 113 by jgs, Mon Feb 28 07:06:33 2005 UTC
# Line 71  void  Finley_RawScaledSystemMatrixVector Line 71  void  Finley_RawScaledSystemMatrixVector
71      double* out) {      double* out) {
72    
73    maybelong ir,icol,iptr,icb,irb,irow,ic;    maybelong ir,icol,iptr,icb,irb,irow,ic;
74      double reg,reg1,reg2,reg3;
75    #pragma omp barrier    #pragma omp barrier
76    
77    if (ABS(beta)>0.) {    if (ABS(beta)>0.) {
# Line 88  void  Finley_RawScaledSystemMatrixVector Line 89  void  Finley_RawScaledSystemMatrixVector
89      if (A ->col_block_size==1 && A->row_block_size ==1) {      if (A ->col_block_size==1 && A->row_block_size ==1) {
90        switch(A->type) {        switch(A->type) {
91        case CSR:        case CSR:
92          #pragma omp for private(irow,iptr) schedule(static)          #pragma omp for private(irow,iptr,reg) schedule(static)
93      for (irow=0;irow< A->pattern->n_ptr;++irow) {      for (irow=0;irow< A->pattern->n_ptr;++irow) {
94              reg=0.;
95        for (iptr=(A->pattern->ptr[irow])-PTR_OFFSET;iptr<(A->pattern->ptr[irow+1])-PTR_OFFSET; ++iptr) {        for (iptr=(A->pattern->ptr[irow])-PTR_OFFSET;iptr<(A->pattern->ptr[irow+1])-PTR_OFFSET; ++iptr) {
96          out[irow] += alpha * A->val[iptr] * in[A->pattern->index[iptr]-INDEX_OFFSET];            reg += A->val[iptr] * in[A->pattern->index[iptr]-INDEX_OFFSET];
97        }        }
98          out[irow] += alpha * reg;
99      }      }
100      break;      break;
101        case CSC:        case CSC:
# Line 108  void  Finley_RawScaledSystemMatrixVector Line 111  void  Finley_RawScaledSystemMatrixVector
111      Finley_ErrorCode=TYPE_ERROR;      Finley_ErrorCode=TYPE_ERROR;
112      sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");      sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");
113        } /* switch A->type */        } /* switch A->type */
114        } else if (A ->col_block_size==2 && A->row_block_size ==2) {
115          switch(A->type) {
116          case CSR:
117            #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg1,reg2) schedule(static)
118        for (ir=0;ir< A->pattern->n_ptr;ir++) {
119              reg1=0.;
120              reg2=0.;
121          for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
122               ic=2*(A->pattern->index[iptr]-INDEX_OFFSET);
123               reg1 += A->val[iptr*4  ]*in[ic] + A->val[iptr*4+2]*in[1+ic];
124               reg2 += A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic];
125          }
126          out[  2*ir] += alpha * reg1;
127          out[1+2*ir] += alpha * reg2;
128        }
129        break;
130          case CSC:
131            /* TODO: parallelize */
132            #pragma omp single
133        for (ic=0;ic< A->pattern->n_ptr;ic++) {
134          for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
135               ic=2*(A->pattern->index[iptr]-INDEX_OFFSET);
136               out[  2*ir] += alpha * ( A->val[iptr*4  ]*in[ic] + A->val[iptr*4+2]*in[1+ic] );
137               out[1+2*ir] += alpha * ( A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic] );
138          }
139        }
140          default:
141        Finley_ErrorCode=TYPE_ERROR;
142        sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");
143          } /* switch A->type */
144        } else if (A ->col_block_size==3 && A->row_block_size ==3) {
145          switch(A->type) {
146          case CSR:
147            #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg1,reg2,reg3) schedule(static)
148        for (ir=0;ir< A->pattern->n_ptr;ir++) {
149              reg1=0.;
150              reg2=0.;
151              reg3=0.;
152          for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
153               ic=3*(A->pattern->index[iptr]-INDEX_OFFSET);
154               reg1 += A->val[iptr*9  ]*in[ic] + A->val[iptr*9+3]*in[1+ic] + A->val[iptr*9+6]*in[2+ic];
155               reg2 += A->val[iptr*9+1]*in[ic] + A->val[iptr*9+4]*in[1+ic] + A->val[iptr*9+7]*in[2+ic];
156               reg3 += A->val[iptr*9+2]*in[ic] + A->val[iptr*9+5]*in[1+ic] + A->val[iptr*9+8]*in[2+ic];
157          }
158          out[  3*ir] += alpha * reg1;
159          out[1+3*ir] += alpha * reg2;
160          out[2+3*ir] += alpha * reg3;
161        }
162        break;
163          case CSC:
164            /* TODO: parallelize */
165            #pragma omp single
166        for (ic=0;ic< A->pattern->n_ptr;ic++) {
167          for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
168              ir=3*(A->pattern->index[iptr]-INDEX_OFFSET);
169                  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] );
170              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] );
171              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] );
172          }
173        }
174          default:
175        Finley_ErrorCode=TYPE_ERROR;
176        sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");
177          } /* switch A->type */
178      } else {      } else {
179        switch(A->type) {        switch(A->type) {
180        case CSR:        case CSR:
181          #pragma omp for private(ir,iptr,irb,icb,irow,icol) schedule(static)          #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg) schedule(static)
182      for (ir=0;ir< A->pattern->n_ptr;ir++) {      for (ir=0;ir< A->pattern->n_ptr;ir++) {
183        for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {        for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
184          for (irb=0;irb< A->row_block_size;irb++) {          for (irb=0;irb< A->row_block_size;irb++) {
185            irow=irb+A->row_block_size*ir;            irow=irb+A->row_block_size*ir;
186                  reg=0.;
187            for (icb=0;icb< A->col_block_size;icb++) {            for (icb=0;icb< A->col_block_size;icb++) {
188          icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);          icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
189          out[irow] += alpha * 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];
190            }            }
191              out[irow] += alpha * reg;
192          }          }
193        }        }
194      }      }
# Line 146  void  Finley_RawScaledSystemMatrixVector Line 215  void  Finley_RawScaledSystemMatrixVector
215    }    }
216    return;    return;
217  }  }
   
 /*  
  * $Log$  
  * Revision 1.4  2004/12/15 07:08:33  jgs  
  * *** empty log message ***  
  *  
  *  
  *  
  */  

Legend:
Removed from v.112  
changed lines
  Added in v.113

  ViewVC Help
Powered by ViewVC 1.1.26