/[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

trunk/esys2/paso/src/SystemMatrix_nullifyRowsAndCols.c revision 150 by jgs, Thu Sep 15 03:44:45 2005 UTC trunk/paso/src/SystemMatrix_nullifyRowsAndCols.c revision 631 by dhawcroft, Thu Mar 23 04:27:32 2006 UTC
# Line 1  Line 1 
1  /* $Id$ */  /* $Id$ */
2    
3    /*
4    ********************************************************************************
5    *               Copyright 2006 by ACcESS MNRF                                *
6    *                                                                              *
7    *                 http://www.access.edu.au                                     *
8    *           Primary Business: Queensland, Australia                            *
9    *     Licensed under the Open Software License version 3.0             *
10    *        http://www.opensource.org/licenses/osl-3.0.php                        *
11    ********************************************************************************
12    */
13    
14  /**************************************************************/  /**************************************************************/
15    
16  /* Paso: SystemMatrix                                       */  /* Paso: SystemMatrix                                       */
# Line 24  Line 35 
35    
36  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) {
37    
38      index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
39    index_t ir,icol,iptr,icb,irb,irow,ic,l;    index_t ir,icol,iptr,icb,irb,irow,ic,l;
40    
41    if (A ->col_block_size==1 && A ->row_block_size ==1) {    if (A ->col_block_size==1 && A ->row_block_size ==1) {
42      switch(A->type) {      if (A->type & MATRIX_FORMAT_CSC) {
     case CSR:  
43        #pragma omp parallel for private(irow, iptr,icol) schedule(static)        #pragma omp parallel for private(irow, iptr,icol) schedule(static)
44        for (irow=0;irow< A->pattern->n_ptr;irow++) {        for (icol=0;icol< A->pattern->n_ptr;icol++) {
45      /* 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++) {
46      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;  
47        if (mask_col[icol]>0. || mask_row[irow]>0. ) {        if (mask_col[icol]>0. || mask_row[irow]>0. ) {
48              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
49            A->val[iptr]=main_diagonal_value;            A->val[iptr]=main_diagonal_value;
# Line 43  void Paso_SystemMatrix_nullifyRowsAndCol Line 53  void Paso_SystemMatrix_nullifyRowsAndCol
53        }        }
54      }      }
55        }        }
56        break;      } else {
     case CSC:  
57        #pragma omp parallel for private(irow, iptr,icol) schedule(static)        #pragma omp parallel for private(irow, iptr,icol) schedule(static)
58        for (icol=0;icol< A->pattern->n_ptr;icol++) {        for (irow=0;irow< A->pattern->n_ptr;irow++) {
59      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 */
60        irow=A->pattern->index[iptr]-INDEX_OFFSET;      for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
61          icol=A->pattern->index[iptr]-index_offset;
62        if (mask_col[icol]>0. || mask_row[irow]>0. ) {        if (mask_col[icol]>0. || mask_row[irow]>0. ) {
63              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
64            A->val[iptr]=main_diagonal_value;            A->val[iptr]=main_diagonal_value;
# Line 58  void Paso_SystemMatrix_nullifyRowsAndCol Line 68  void Paso_SystemMatrix_nullifyRowsAndCol
68        }        }
69      }      }
70        }        }
71        break;      }
     default:  
       Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");  
     } /* switch A->type */  
72    } else {    } else {
73      switch(A->type) {      if (A->type & MATRIX_FORMAT_CSC) {
74      case CSR:        #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
75        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)        for (ic=0;ic< A->pattern->n_ptr;ic++) {
76        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++) {  
77        for (irb=0;irb< A->row_block_size;irb++) {        for (irb=0;irb< A->row_block_size;irb++) {
78          irow=irb+A->row_block_size*ir;          irow=irb+A->row_block_size*(A->pattern->index[iptr]-index_offset);
79          for (icb=0;icb< A->col_block_size;icb++) {          for (icb=0;icb< A->col_block_size;icb++) {
80            icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);            icol=icb+A->col_block_size*ic;
81            if (mask_col[icol]>0. || mask_row[irow]>0. ) {            if (mask_col[icol]>0. || mask_row[irow]>0. ) {
82                  l=iptr*A->block_size+irb+A->row_block_size*icb;                  l=iptr*A->block_size+irb+A->row_block_size*icb;
83          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 90  void Paso_SystemMatrix_nullifyRowsAndCol
90        }        }
91      }      }
92        }        }
93        break;      } else {
94      case CSC:        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
95        #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)        for (ir=0;ir< A->pattern->n_ptr;ir++) {
96        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++) {  
97        for (irb=0;irb< A->row_block_size;irb++) {        for (irb=0;irb< A->row_block_size;irb++) {
98          irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);          irow=irb+A->row_block_size*ir;
99          for (icb=0;icb< A->col_block_size;icb++) {          for (icb=0;icb< A->col_block_size;icb++) {
100            icol=icb+A->col_block_size*ic;            icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
101            if (mask_col[icol]>0. || mask_row[irow]>0. ) {            if (mask_col[icol]>0. || mask_row[irow]>0. ) {
102                  l=iptr*A->block_size+irb+A->row_block_size*icb;                  l=iptr*A->block_size+irb+A->row_block_size*icb;
103          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 110  void Paso_SystemMatrix_nullifyRowsAndCol
110        }        }
111      }      }
112        }        }
113        break;     }
     default:  
       Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");  
     } /* switch A->type */  
114    }    }
115    return;    return;
116  }  }

Legend:
Removed from v.150  
changed lines
  Added in v.631

  ViewVC Help
Powered by ViewVC 1.1.26