/[escript]/trunk/paso/src/SystemMatrix_nullifyRowsAndCols.c
ViewVC logotype

Diff of /trunk/paso/src/SystemMatrix_nullifyRowsAndCols.c

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

revision 155 by jgs, Wed Nov 9 02:02:19 2005 UTC revision 415 by gross, Wed Jan 4 05:37:33 2006 UTC
# Line 24  Line 24 
24    
25  void Paso_SystemMatrix_nullifyRowsAndCols(Paso_SystemMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {  void Paso_SystemMatrix_nullifyRowsAndCols(Paso_SystemMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
26    
27      index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
28    index_t ir,icol,iptr,icb,irb,irow,ic,l;    index_t ir,icol,iptr,icb,irb,irow,ic,l;
29    
30    if (A ->col_block_size==1 && A ->row_block_size ==1) {    if (A ->col_block_size==1 && A ->row_block_size ==1) {
31      switch(A->type) {      if (A->type & MATRIX_FORMAT_CSC) {
     case CSR:  
32        #pragma omp parallel for private(irow, iptr,icol) schedule(static)        #pragma omp parallel for private(irow, iptr,icol) schedule(static)
33        for (irow=0;irow< A->pattern->n_ptr;irow++) {        for (icol=0;icol< A->pattern->n_ptr;icol++) {
34      /* TODO: test mask_row here amd not inside every loop */      for (iptr=A->pattern->ptr[icol]-index_offset;iptr<A->pattern->ptr[icol+1]-index_offset; iptr++) {
35      for (iptr=A->pattern->ptr[irow]-PTR_OFFSET;iptr<A->pattern->ptr[irow+1]-PTR_OFFSET; iptr++) {        irow=A->pattern->index[iptr]-index_offset;
       icol=A->pattern->index[iptr]-INDEX_OFFSET;  
36        if (mask_col[icol]>0. || mask_row[irow]>0. ) {        if (mask_col[icol]>0. || mask_row[irow]>0. ) {
37              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
38            A->val[iptr]=main_diagonal_value;            A->val[iptr]=main_diagonal_value;
# Line 43  void Paso_SystemMatrix_nullifyRowsAndCol Line 42  void Paso_SystemMatrix_nullifyRowsAndCol
42        }        }
43      }      }
44        }        }
45        break;      } else {
     case CSC:  
46        #pragma omp parallel for private(irow, iptr,icol) schedule(static)        #pragma omp parallel for private(irow, iptr,icol) schedule(static)
47        for (icol=0;icol< A->pattern->n_ptr;icol++) {        for (irow=0;irow< A->pattern->n_ptr;irow++) {
48      for (iptr=A->pattern->ptr[icol]-PTR_OFFSET;iptr<A->pattern->ptr[icol+1]-PTR_OFFSET; iptr++) {      /* TODO: test mask_row here amd not inside every loop */
49        irow=A->pattern->index[iptr]-INDEX_OFFSET;      for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
50          icol=A->pattern->index[iptr]-index_offset;
51        if (mask_col[icol]>0. || mask_row[irow]>0. ) {        if (mask_col[icol]>0. || mask_row[irow]>0. ) {
52              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
53            A->val[iptr]=main_diagonal_value;            A->val[iptr]=main_diagonal_value;
# Line 58  void Paso_SystemMatrix_nullifyRowsAndCol Line 57  void Paso_SystemMatrix_nullifyRowsAndCol
57        }        }
58      }      }
59        }        }
60        break;      }
     default:  
       Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");  
     } /* switch A->type */  
61    } else {    } else {
62      switch(A->type) {      if (A->type & MATRIX_FORMAT_CSC) {
63      case CSR:        #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
64        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)        for (ic=0;ic< A->pattern->n_ptr;ic++) {
65        for (ir=0;ir< A->pattern->n_ptr;ir++) {      for (iptr=A->pattern->ptr[ic]-index_offset;iptr<A->pattern->ptr[ic+1]-index_offset; iptr++) {
     for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {  
66        for (irb=0;irb< A->row_block_size;irb++) {        for (irb=0;irb< A->row_block_size;irb++) {
67          irow=irb+A->row_block_size*ir;          irow=irb+A->row_block_size*(A->pattern->index[iptr]-index_offset);
68          for (icb=0;icb< A->col_block_size;icb++) {          for (icb=0;icb< A->col_block_size;icb++) {
69            icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);            icol=icb+A->col_block_size*ic;
70            if (mask_col[icol]>0. || mask_row[irow]>0. ) {            if (mask_col[icol]>0. || mask_row[irow]>0. ) {
71                  l=iptr*A->block_size+irb+A->row_block_size*icb;                  l=iptr*A->block_size+irb+A->row_block_size*icb;
72          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
# Line 84  void Paso_SystemMatrix_nullifyRowsAndCol Line 79  void Paso_SystemMatrix_nullifyRowsAndCol
79        }        }
80      }      }
81        }        }
82        break;      } else {
83      case CSC:        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
84        #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)        for (ir=0;ir< A->pattern->n_ptr;ir++) {
85        for (ic=0;ic< A->pattern->n_ptr;ic++) {      for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
     for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {  
86        for (irb=0;irb< A->row_block_size;irb++) {        for (irb=0;irb< A->row_block_size;irb++) {
87          irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);          irow=irb+A->row_block_size*ir;
88          for (icb=0;icb< A->col_block_size;icb++) {          for (icb=0;icb< A->col_block_size;icb++) {
89            icol=icb+A->col_block_size*ic;            icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
90            if (mask_col[icol]>0. || mask_row[irow]>0. ) {            if (mask_col[icol]>0. || mask_row[irow]>0. ) {
91                  l=iptr*A->block_size+irb+A->row_block_size*icb;                  l=iptr*A->block_size+irb+A->row_block_size*icb;
92          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
# Line 105  void Paso_SystemMatrix_nullifyRowsAndCol Line 99  void Paso_SystemMatrix_nullifyRowsAndCol
99        }        }
100      }      }
101        }        }
102        break;     }
     default:  
       Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");  
     } /* switch A->type */  
103    }    }
104    return;    return;
105  }  }

Legend:
Removed from v.155  
changed lines
  Added in v.415

  ViewVC Help
Powered by ViewVC 1.1.26