/[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 99 by jgs, Tue Dec 14 05:39:33 2004 UTC revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC
# Line 46  void  Finley_ScaledSystemMatrixVector(do Line 46  void  Finley_ScaledSystemMatrixVector(do
46    } else if (getLength(out) != A ->row_block_size * A ->num_rows) {    } else if (getLength(out) != A ->row_block_size * A ->num_rows) {
47      Finley_ErrorCode=TYPE_ERROR;      Finley_ErrorCode=TYPE_ERROR;
48      sprintf(Finley_ErrorMsg,"matrix vector product : lnegth of output Data object and matrix row dimensions don't match.");      sprintf(Finley_ErrorMsg,"matrix vector product : lnegth of output Data object and matrix row dimensions don't match.");
49      } else if (A ->symmetric) {
50        Finley_ErrorCode=SYSTEM_ERROR;
51        sprintf(Finley_ErrorMsg,"matrix-vector product for symmetrix matric has not been implemented yet.");
52    }    }
53    
54    /* delegate to routine */    /* delegate to routine */
# Line 89  void  Finley_RawScaledSystemMatrixVector Line 92  void  Finley_RawScaledSystemMatrixVector
92        switch(A->type) {        switch(A->type) {
93        case CSR:        case CSR:
94          #pragma omp for private(irow,iptr) schedule(static)          #pragma omp for private(irow,iptr) schedule(static)
95      for (irow=0;irow< A->pattern->n_ptr;++irow) {      for (irow=0;irow< A->num_rows;irow++) {
96        for (iptr=(A->pattern->ptr[irow])-PTR_OFFSET;iptr<(A->pattern->ptr[irow+1])-PTR_OFFSET; ++iptr) {        for (iptr=A->ptr[irow]-PTR_OFFSET;iptr<A->ptr[irow+1]-PTR_OFFSET; iptr++) {
97          out[irow] += alpha * A->val[iptr] * in[A->pattern->index[iptr]-INDEX_OFFSET];          out[irow] += alpha * A->val[iptr] * in[A->index[iptr]-INDEX_OFFSET];
98        }        }
99      }      }
100      break;      break;
101        case CSC:        case CSC:
102          /* TODO: parallelize */          /* TODO: parallelize */
103          #pragma omp single          #pragma omp single
104      for (icol=0;icol< A->pattern->n_ptr;++icol) {      for (icol=0;icol< A->num_cols;icol++) {
105        for (iptr=A->pattern->ptr[icol]-PTR_OFFSET;iptr<A->pattern->ptr[icol+1]-PTR_OFFSET; ++iptr) {        for (iptr=A->ptr[icol]-PTR_OFFSET;iptr<A->ptr[icol+1]-PTR_OFFSET; iptr++) {
106          out[A->pattern->index[iptr]-INDEX_OFFSET]+= alpha * A->val[iptr] * in[icol];          out[A->index[iptr]-INDEX_OFFSET]+= alpha * A->val[iptr] * in[icol];
107        }        }
108      }      }
109      break;      break;
110        default:        default:
111      Finley_ErrorCode=TYPE_ERROR;      Finley_ErrorCode=TYPE_ERROR;
112      sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");      sprintf(Finley_ErrorMsg,"Unknown matrix type.");
113        } /* switch A->type */        } /* switch A->type */
114      } else {      } else {
115        switch(A->type) {        switch(A->type) {
116        case CSR:        case CSR:
117          #pragma omp for private(ir,iptr,irb,icb,irow,icol) schedule(static)          #pragma omp for private(ir,iptr,irb,icb,irow,icol) schedule(static)
118      for (ir=0;ir< A->pattern->n_ptr;ir++) {      for (ir=0;ir< A->num_rows;ir++) {
119        for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {        for (iptr=A->ptr[ir]-PTR_OFFSET;iptr<A->ptr[ir+1]-PTR_OFFSET; iptr++) {
120          for (irb=0;irb< A->row_block_size;irb++) {          for (irb=0;irb< A->row_block_size;irb++) {
121            irow=irb+A->row_block_size*ir;            irow=irb+A->row_block_size*ir;
122            for (icb=0;icb< A->col_block_size;icb++) {            for (icb=0;icb< A->col_block_size;icb++) {
123          icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);          icol=icb+A->col_block_size*(A->index[iptr]-INDEX_OFFSET);
124          out[irow] += alpha * A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];          out[irow] += alpha * A->val[iptr*A->row_block_size*A->col_block_size+icb+A->col_block_size*irb] * in[icol];
125            }            }
126          }          }
127        }        }
# Line 127  void  Finley_RawScaledSystemMatrixVector Line 130  void  Finley_RawScaledSystemMatrixVector
130        case CSC:        case CSC:
131          /* TODO: parallelize */          /* TODO: parallelize */
132          #pragma omp single          #pragma omp single
133      for (ic=0;ic< A->pattern->n_ptr;ic++) {      for (ic=0;ic< A->num_cols;ic++) {
134        for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {        for (iptr=A->ptr[ic]-PTR_OFFSET;iptr<A->ptr[ic+1]-PTR_OFFSET; iptr++) {
135          for (irb=0;irb< A->row_block_size;irb++) {          for (irb=0;irb< A->row_block_size;irb++) {
136            irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);            irow=irb+A->row_block_size*(A->index[iptr]-INDEX_OFFSET);
137            for (icb=0;icb< A->col_block_size;icb++) {            for (icb=0;icb< A->col_block_size;icb++) {
138          icol=icb+A->col_block_size*ic;          icol=icb+A->col_block_size*ic;
139          out[irow] += alpha * A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];          out[irow] += alpha * A->val[iptr*A->row_block_size*A->col_block_size+icb+A->col_block_size*irb] * in[icol];
140            }            }
141          }          }
142        }        }
143      }      }
144        default:        default:
145      Finley_ErrorCode=TYPE_ERROR;      Finley_ErrorCode=TYPE_ERROR;
146      sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");      sprintf(Finley_ErrorMsg,"Unknown matrix type.");
147        } /* switch A->type */        } /* switch A->type */
148      }      }
149    }    }
# Line 149  void  Finley_RawScaledSystemMatrixVector Line 152  void  Finley_RawScaledSystemMatrixVector
152    
153  /*  /*
154   * $Log$   * $Log$
155   * Revision 1.2  2004/12/14 05:39:30  jgs   * Revision 1.3  2004/12/15 03:48:46  jgs
156   * *** empty log message ***   * *** empty log message ***
157   *   *
  * Revision 1.1.1.1.2.2  2004/12/07 10:12:05  gross  
  * GMRES added  
  *  
  * Revision 1.1.1.1.2.1  2004/11/12 06:58:19  gross  
  * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry  
  *  
158   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs
159   * initial import of project esys2   * initial import of project esys2
160   *   *

Legend:
Removed from v.99  
changed lines
  Added in v.100

  ViewVC Help
Powered by ViewVC 1.1.26