/[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 969 by ksteube, Tue Feb 13 23:02:23 2007 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 if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
57      case CSC:        /*
58            Plan: Collect non-zeros on each CPU,
59        use MPI to concatenate,
60        sort with qsort,
61        and then use binary search to check if a row is non-zero
62          */
63          fprintf(stderr, "SystemMatrix_nullifyRowsAndCols: need to implement MATRIX_FORMAT_TRILINOS_CRS (and for block sizes != 1)\n");
64          exit(1);
65        } else {
66        #pragma omp parallel for private(irow, iptr,icol) schedule(static)        #pragma omp parallel for private(irow, iptr,icol) schedule(static)
67        for (icol=0;icol< A->pattern->n_ptr;icol++) {        for (irow=0;irow< A->pattern->n_ptr;irow++) {
68      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 */
69        irow=A->pattern->index[iptr]-INDEX_OFFSET;      for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
70          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 58  void Paso_SystemMatrix_nullifyRowsAndCol Line 77  void Paso_SystemMatrix_nullifyRowsAndCol
77        }        }
78      }      }
79        }        }
80        break;      }
     default:  
       Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");  
     } /* switch A->type */  
81    } else {    } else {
82      switch(A->type) {      if (A->type & MATRIX_FORMAT_CSC) {
83      case CSR:        #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
84        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)        for (ic=0;ic< A->pattern->n_ptr;ic++) {
85        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++) {  
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*ir;          irow=irb+A->row_block_size*(A->pattern->index[iptr]-index_offset);
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*(A->pattern->index[iptr]-INDEX_OFFSET);            icol=icb+A->col_block_size*ic;
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 84  void Paso_SystemMatrix_nullifyRowsAndCol Line 99  void Paso_SystemMatrix_nullifyRowsAndCol
99        }        }
100      }      }
101        }        }
102        break;      } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
103      case CSC:        fprintf(stderr, "SystemMatrix_nullifyRowsAndCols 2: need to implement MATRIX_FORMAT_TRILINOS_CRS\n");
104        #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)        exit(1);
105        for (ic=0;ic< A->pattern->n_ptr;ic++) {      } else {
106      for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
107          for (ir=0;ir< A->pattern->n_ptr;ir++) {
108        for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
109        for (irb=0;irb< A->row_block_size;irb++) {        for (irb=0;irb< A->row_block_size;irb++) {
110          irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);          irow=irb+A->row_block_size*ir;
111          for (icb=0;icb< A->col_block_size;icb++) {          for (icb=0;icb< A->col_block_size;icb++) {
112            icol=icb+A->col_block_size*ic;            icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
113            if (mask_col[icol]>0. || mask_row[irow]>0. ) {            if (mask_col[icol]>0. || mask_row[irow]>0. ) {
114                  l=iptr*A->block_size+irb+A->row_block_size*icb;                  l=iptr*A->block_size+irb+A->row_block_size*icb;
115          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 122  void Paso_SystemMatrix_nullifyRowsAndCol
122        }        }
123      }      }
124        }        }
125        break;     }
     default:  
       Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");  
     } /* switch A->type */  
126    }    }
127    return;    return;
128  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26