/[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 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC
# Line 1  Line 1 
1  /* $Id$ */  
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2008 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
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    
17  /* Paso: SystemMatrix                                       */  /* Paso: SystemMatrix                                       */
18    
19  /*  nullify rows and columns in the matrix                    */  /*  nullify rows and columns in the matrix                    */
20    /*
21  /*  the rows and columns are marked by positive values in     */  /*  the rows and columns are marked by positive values in     */
22  /*  mask_row and mask_col. Values on the main diagonal        */  /*  mask_row and mask_col. Values on the main diagonal        */
23  /*  which are marked to set to zero by both mask_row and      */  /*  which are marked to set to zero by both mask_row and      */
# Line 24  Line 36 
36    
37  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) {
38    
39    index_t ir,icol,iptr,icb,irb,irow,ic,l;    double *remote_values=NULL;
40      Paso_MPIInfo *mpi_info=A->mpi_info;
41    if (A ->col_block_size==1 && A ->row_block_size ==1) {    if (mpi_info->size>1) {
42      switch(A->type) {       if (A ->col_block_size==1 && A ->row_block_size ==1) {
43      case CSR:         if (A->type & MATRIX_FORMAT_CSC) {
44        #pragma omp parallel for private(irow, iptr,icol) schedule(static)             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: CSC is not supported by MPI.");
45        for (irow=0;irow< A->pattern->n_ptr;irow++) {             return;
46      /* TODO: test mask_row here amd not inside every loop */         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
47      for (iptr=A->pattern->ptr[irow]-PTR_OFFSET;iptr<A->pattern->ptr[irow+1]-PTR_OFFSET; iptr++) {             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported with MPI.");
48        icol=A->pattern->index[iptr]-INDEX_OFFSET;             return;
49        if (mask_col[icol]>0. || mask_row[irow]>0. ) {         } else {
50              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {           if (Paso_noError()) {
51            A->val[iptr]=main_diagonal_value;              Paso_SystemMatrix_startColCollect(A,mask_col);
52              } else {              Paso_SystemMatrix_startRowCollect(A,mask_row);
53            A->val[iptr]=0;              Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->mainBlock,mask_row,mask_col,main_diagonal_value);
54              }              remote_values=Paso_SystemMatrix_finishColCollect(A);
55        }              Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->col_coupleBlock,mask_row,remote_values,0.);
56      }              remote_values=Paso_SystemMatrix_finishRowCollect(A);
57        }              Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->row_coupleBlock,remote_values,mask_col,0.);
58        break;           }
59      case CSC:         }
60        #pragma omp parallel for private(irow, iptr,icol) schedule(static)       } else {
61        for (icol=0;icol< A->pattern->n_ptr;icol++) {         if (A->type & MATRIX_FORMAT_CSC) {
62      for (iptr=A->pattern->ptr[icol]-PTR_OFFSET;iptr<A->pattern->ptr[icol+1]-PTR_OFFSET; iptr++) {             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: CSC is not supported by MPI.");
63        irow=A->pattern->index[iptr]-INDEX_OFFSET;             return;
64        if (mask_col[icol]>0. || mask_row[irow]>0. ) {         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
65              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported with MPI.");
66            A->val[iptr]=main_diagonal_value;             return;
67              } else {         } else {
68            A->val[iptr]=0;           if (Paso_noError()) {
69              }              Paso_SystemMatrix_startColCollect(A,mask_col) ;
70        }              Paso_SystemMatrix_startRowCollect(A,mask_row) ;
71      }              Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->mainBlock,mask_row,mask_col,main_diagonal_value);
72        }              remote_values=Paso_SystemMatrix_finishColCollect(A);
73        break;              Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->col_coupleBlock,mask_row,remote_values,0.);
74      default:              remote_values=Paso_SystemMatrix_finishRowCollect(A);
75        Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");              Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->row_coupleBlock,remote_values,mask_col,0.);
76      } /* switch A->type */           }
77    } else {         }
78      switch(A->type) {       }
79      case CSR:    } else {
80        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)       if (A ->col_block_size==1 && A ->row_block_size ==1) {
81        for (ir=0;ir< A->pattern->n_ptr;ir++) {         if (A->type & MATRIX_FORMAT_CSC) {
82      for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {             Paso_SparseMatrix_nullifyRowsAndCols_CSC_BLK1(A->mainBlock,mask_row,mask_col,main_diagonal_value);
83        for (irb=0;irb< A->row_block_size;irb++) {         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
84          irow=irb+A->row_block_size*ir;             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported.");
85          for (icb=0;icb< A->col_block_size;icb++) {         } else {
86            icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);             Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->mainBlock,mask_row,mask_col,main_diagonal_value);
87            if (mask_col[icol]>0. || mask_row[irow]>0. ) {         }
88                  l=iptr*A->block_size+irb+A->row_block_size*icb;       } else {
89          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {         if (A->type & MATRIX_FORMAT_CSC) {
90            A->val[l]=main_diagonal_value;             Paso_SparseMatrix_nullifyRowsAndCols_CSC(A->mainBlock,mask_row,mask_col,main_diagonal_value);
91          } else {         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
92            A->val[l]=0;             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported with MPI.");
93          }         } else {
94            }             Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->mainBlock,mask_row,mask_col,main_diagonal_value);
95          }         }
96        }       }
     }  
       }  
       break;  
     case CSC:  
       #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)  
       for (ic=0;ic< A->pattern->n_ptr;ic++) {  
     for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {  
       for (irb=0;irb< A->row_block_size;irb++) {  
         irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);  
         for (icb=0;icb< A->col_block_size;icb++) {  
           icol=icb+A->col_block_size*ic;  
           if (mask_col[icol]>0. || mask_row[irow]>0. ) {  
                 l=iptr*A->block_size+irb+A->row_block_size*icb;  
         if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {  
           A->val[l]=main_diagonal_value;  
         } else {  
           A->val[l]=0;  
         }  
           }  
         }  
       }  
     }  
       }  
       break;  
     default:  
       Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");  
     } /* switch A->type */  
97    }    }
98    return;    return;
99  }  }
   
 /*  
  * $Log$  
  * Revision 1.2  2005/09/15 03:44:39  jgs  
  * Merge of development branch dev-02 back to main trunk on 2005-09-15  
  *  
  * Revision 1.1.2.1  2005/09/05 06:29:48  gross  
  * These files have been extracted from finley to define a stand alone libray for iterative  
  * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but  
  * has not been tested yet.  
  *  
  *  
  */  

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

  ViewVC Help
Powered by ViewVC 1.1.26