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

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

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

revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC
# Line 58  void  Finley_SystemMatrix_nullifyRowsAnd Line 58  void  Finley_SystemMatrix_nullifyRowsAnd
58    
59  void Finley_SystemMatrixNullify(Finley_SystemMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {  void Finley_SystemMatrixNullify(Finley_SystemMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
60    
61    maybelong ir,icol,iptr,icb,irb,irow,ic,l;    maybelong ir,icol,iptr,icb,irb,irow,ic;
62    
63    if (A ->col_block_size==1 && A ->row_block_size ==1) {    if (A ->col_block_size==1 && A ->row_block_size ==1) {
64      switch(A->type) {      switch(A->type) {
65      case CSR:      case CSR:
66        #pragma omp parallel for private(irow, iptr,icol) schedule(static)        #pragma omp parallel for private(irow, iptr,icol) schedule(static)
67        for (irow=0;irow< A->pattern->n_ptr;irow++) {        for (irow=0;irow< A->num_rows;irow++) {
68      /* TODO: test mask_row here amd not inside every loop */      /* TODO: test mask_row here amd not inside every loop */
69      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++) {
70        icol=A->pattern->index[iptr]-INDEX_OFFSET;        icol=A->index[iptr]-INDEX_OFFSET;
71        if (mask_col[icol]>0. || mask_row[irow]>0. ) {        if (mask_col[icol]>0. || mask_row[irow]>0. ) {
72              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
73            A->val[iptr]=main_diagonal_value;            A->val[iptr]=main_diagonal_value;
# Line 80  void Finley_SystemMatrixNullify(Finley_S Line 80  void Finley_SystemMatrixNullify(Finley_S
80        break;        break;
81      case CSC:      case CSC:
82        #pragma omp parallel for private(irow, iptr,icol) schedule(static)        #pragma omp parallel for private(irow, iptr,icol) schedule(static)
83        for (icol=0;icol< A->pattern->n_ptr;icol++) {        for (icol=0;icol< A->num_cols;icol++) {
84      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++) {
85        irow=A->pattern->index[iptr]-INDEX_OFFSET;        irow=A->index[iptr]-INDEX_OFFSET;
86        if (mask_col[icol]>0. || mask_row[irow]>0. ) {        if (mask_col[icol]>0. || mask_row[irow]>0. ) {
87              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
88            A->val[iptr]=main_diagonal_value;            A->val[iptr]=main_diagonal_value;
# Line 100  void Finley_SystemMatrixNullify(Finley_S Line 100  void Finley_SystemMatrixNullify(Finley_S
100    } else {    } else {
101      switch(A->type) {      switch(A->type) {
102      case CSR:      case CSR:
103        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)        #pragma omp parallel for private(irow, iptr,icol,ir,irb,icb) schedule(static)
104        for (ir=0;ir< A->pattern->n_ptr;ir++) {        for (ir=0;ir< A->num_rows;ir++) {
105      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++) {
106        for (irb=0;irb< A->row_block_size;irb++) {        for (irb=0;irb< A->row_block_size;irb++) {
107          irow=irb+A->row_block_size*ir;          irow=irb+A->row_block_size*ir;
108          for (icb=0;icb< A->col_block_size;icb++) {          for (icb=0;icb< A->col_block_size;icb++) {
109            icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);            icol=icb+A->col_block_size*(A->index[iptr]-INDEX_OFFSET);
110            if (mask_col[icol]>0. || mask_row[irow]>0. ) {            if (mask_col[icol]>0. || mask_row[irow]>0. ) {
                 l=iptr*A->block_size+irb+A->row_block_size*icb;  
111          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
112            A->val[l]=main_diagonal_value;            A->val[iptr]=main_diagonal_value;
113          } else {          } else {
114            A->val[l]=0;            A->val[iptr]=0;
115          }          }
116            }            }
117          }          }
# Line 121  void Finley_SystemMatrixNullify(Finley_S Line 120  void Finley_SystemMatrixNullify(Finley_S
120        }        }
121        break;        break;
122      case CSC:      case CSC:
123        #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)        #pragma omp parallel for private(irow, iptr,icol,ic,irb,icb) schedule(static)
124        for (ic=0;ic< A->pattern->n_ptr;ic++) {        for (ic=0;ic< A->num_cols;ic++) {
125      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++) {
126        for (irb=0;irb< A->row_block_size;irb++) {        for (irb=0;irb< A->row_block_size;irb++) {
127          irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);          irow=irb+A->row_block_size*(A->index[iptr]-INDEX_OFFSET);
128          for (icb=0;icb< A->col_block_size;icb++) {          for (icb=0;icb< A->col_block_size;icb++) {
129            icol=icb+A->col_block_size*ic;            icol=icb+A->col_block_size*ic;
130            if (mask_col[icol]>0. || mask_row[irow]>0. ) {            if (mask_col[icol]>0. || mask_row[irow]>0. ) {
                 l=iptr*A->block_size+irb+A->row_block_size*icb;  
131          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
132            A->val[l]=main_diagonal_value;            A->val[iptr]=main_diagonal_value;
133          } else {          } else {
134            A->val[l]=0;            A->val[iptr]=0;
135          }          }
136            }            }
137          }          }
# Line 151  void Finley_SystemMatrixNullify(Finley_S Line 149  void Finley_SystemMatrixNullify(Finley_S
149    
150  /*  /*
151   * $Log$   * $Log$
152   * Revision 1.2  2004/12/14 05:39:31  jgs   * Revision 1.3  2004/12/15 03:48:46  jgs
153   * *** empty log message ***   * *** empty log message ***
154   *   *
  * 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  
  *  
155   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs
156   * initial import of project esys2   * initial import of project esys2
157   *   *

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

  ViewVC Help
Powered by ViewVC 1.1.26