/[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 100 by jgs, Wed Dec 15 03:48:48 2004 UTC revision 102 by jgs, Wed Dec 15 07:08:39 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;    maybelong ir,icol,iptr,icb,irb,irow,ic,l;
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->num_rows;irow++) {        for (irow=0;irow< A->pattern->n_ptr;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->ptr[irow]-PTR_OFFSET;iptr<A->ptr[irow+1]-PTR_OFFSET; iptr++) {      for (iptr=A->pattern->ptr[irow]-PTR_OFFSET;iptr<A->pattern->ptr[irow+1]-PTR_OFFSET; iptr++) {
70        icol=A->index[iptr]-INDEX_OFFSET;        icol=A->pattern->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->num_cols;icol++) {        for (icol=0;icol< A->pattern->n_ptr;icol++) {
84      for (iptr=A->ptr[icol]-PTR_OFFSET;iptr<A->ptr[icol+1]-PTR_OFFSET; iptr++) {      for (iptr=A->pattern->ptr[icol]-PTR_OFFSET;iptr<A->pattern->ptr[icol+1]-PTR_OFFSET; iptr++) {
85        irow=A->index[iptr]-INDEX_OFFSET;        irow=A->pattern->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(irow, iptr,icol,ir,irb,icb) schedule(static)        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
104        for (ir=0;ir< A->num_rows;ir++) {        for (ir=0;ir< A->pattern->n_ptr;ir++) {
105      for (iptr=A->ptr[ir]-PTR_OFFSET;iptr<A->ptr[ir+1]-PTR_OFFSET; iptr++) {      for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->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->index[iptr]-INDEX_OFFSET);            icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
110            if (mask_col[icol]>0. || mask_row[irow]>0. ) {            if (mask_col[icol]>0. || mask_row[irow]>0. ) {
111                    l=iptr*A->block_size+irb+A->row_block_size*icb;
112          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
113            A->val[iptr]=main_diagonal_value;            A->val[l]=main_diagonal_value;
114          } else {          } else {
115            A->val[iptr]=0;            A->val[l]=0;
116          }          }
117            }            }
118          }          }
# Line 120  void Finley_SystemMatrixNullify(Finley_S Line 121  void Finley_SystemMatrixNullify(Finley_S
121        }        }
122        break;        break;
123      case CSC:      case CSC:
124        #pragma omp parallel for private(irow, iptr,icol,ic,irb,icb) schedule(static)        #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
125        for (ic=0;ic< A->num_cols;ic++) {        for (ic=0;ic< A->pattern->n_ptr;ic++) {
126      for (iptr=A->ptr[ic]-PTR_OFFSET;iptr<A->ptr[ic+1]-PTR_OFFSET; iptr++) {      for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
127        for (irb=0;irb< A->row_block_size;irb++) {        for (irb=0;irb< A->row_block_size;irb++) {
128          irow=irb+A->row_block_size*(A->index[iptr]-INDEX_OFFSET);          irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
129          for (icb=0;icb< A->col_block_size;icb++) {          for (icb=0;icb< A->col_block_size;icb++) {
130            icol=icb+A->col_block_size*ic;            icol=icb+A->col_block_size*ic;
131            if (mask_col[icol]>0. || mask_row[irow]>0. ) {            if (mask_col[icol]>0. || mask_row[irow]>0. ) {
132                    l=iptr*A->block_size+irb+A->row_block_size*icb;
133          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
134            A->val[iptr]=main_diagonal_value;            A->val[l]=main_diagonal_value;
135          } else {          } else {
136            A->val[iptr]=0;            A->val[l]=0;
137          }          }
138            }            }
139          }          }
# Line 149  void Finley_SystemMatrixNullify(Finley_S Line 151  void Finley_SystemMatrixNullify(Finley_S
151    
152  /*  /*
153   * $Log$   * $Log$
154   * Revision 1.3  2004/12/15 03:48:46  jgs   * Revision 1.4  2004/12/15 07:08:33  jgs
155   * *** empty log message ***   * *** empty log message ***
156   *   *
  * Revision 1.1.1.1  2004/10/26 06:53:57  jgs  
  * initial import of project esys2  
  *  
  * Revision 1.1  2004/07/02 04:21:13  gross  
  * Finley C code has been included  
  *  
157   *   *
158   */   */

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

  ViewVC Help
Powered by ViewVC 1.1.26