/[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 971 by ksteube, Wed Feb 14 04:40:49 2007 UTC revision 1563 by gross, Wed May 21 15:27:34 2008 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /* $Id$ */
3    
4  /*  /*******************************************************
5  ********************************************************************************   *
6  *               Copyright   2006 by ACcESS MNRF                                *   *           Copyright 2003-2007 by ACceSS MNRF
7  *                                                                              *   *       Copyright 2007 by University of Queensland
8  *                 http://www.access.edu.au                                     *   *
9  *           Primary Business: Queensland, Australia                            *   *                http://esscc.uq.edu.au
10  *     Licensed under the Open Software License version 3.0             *   *        Primary Business: Queensland, Australia
11  *        http://www.opensource.org/licenses/osl-3.0.php                        *   *  Licensed under the Open Software License version 3.0
12  ********************************************************************************   *     http://www.opensource.org/licenses/osl-3.0.php
13  */   *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
18  /* Paso: SystemMatrix                                       */  /* Paso: SystemMatrix                                       */
19    
20  /*  nullify rows and columns in the matrix                    */  /*  nullify rows and columns in the matrix                    */
21    /*
22  /*  the rows and columns are marked by positive values in     */  /*  the rows and columns are marked by positive values in     */
23  /*  mask_row and mask_col. Values on the main diagonal        */  /*  mask_row and mask_col. Values on the main diagonal        */
24  /*  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 35  Line 37 
37    
38  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) {
39    
40    index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);    double *remote_values=NULL;
41    index_t ir,icol,iptr,icb,irb,irow,ic,l;    Paso_MPIInfo *mpi_info=A->mpi_info;
42      if (mpi_info->size>1) {
43    if (A ->col_block_size==1 && A ->row_block_size ==1) {       if (A ->col_block_size==1 && A ->row_block_size ==1) {
44      if (A->type & MATRIX_FORMAT_CSC) {         if (A->type & MATRIX_FORMAT_CSC) {
45        #pragma omp parallel for private(irow, iptr,icol) schedule(static)             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: CSC is not supported by MPI.");
46        for (icol=0;icol< A->pattern->n_ptr;icol++) {             return;
47      for (iptr=A->pattern->ptr[icol]-index_offset;iptr<A->pattern->ptr[icol+1]-index_offset; iptr++) {         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
48        irow=A->pattern->index[iptr]-index_offset;             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported with MPI.");
49        if (mask_col[icol]>0. || mask_row[irow]>0. ) {             return;
50              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {         } else {
51            A->val[iptr]=main_diagonal_value;           if (Paso_noError()) {
52              } else {              Paso_SystemMatrix_startColCollect(A,mask_col);
53            A->val[iptr]=0;              Paso_SystemMatrix_startRowCollect(A,mask_row);
54              }              Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->mainBlock,mask_row,mask_col,main_diagonal_value);
55        }              remote_values=Paso_SystemMatrix_finishColCollect(A);
56      }              Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->col_coupleBlock,mask_row,remote_values,0.);
57        }              remote_values=Paso_SystemMatrix_finishRowCollect(A);
58      } else {              Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->row_coupleBlock,remote_values,mask_col,0.);
59        #pragma omp parallel for private(irow, iptr,icol) schedule(static)           }
60        for (irow=0;irow< A->pattern->n_ptr;irow++) {         }
61      /* TODO: test mask_row here amd not inside every loop */       } else {
62      for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {         if (A->type & MATRIX_FORMAT_CSC) {
63        icol=A->pattern->index[iptr]-index_offset;             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: CSC is not supported by MPI.");
64        if (mask_col[icol]>0. || mask_row[irow]>0. ) {             return;
65              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
66            A->val[iptr]=main_diagonal_value;             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported with MPI.");
67              } else {             return;
68            A->val[iptr]=0;         } else {
69              }           if (Paso_noError()) {
70        }              Paso_SystemMatrix_startColCollect(A,mask_col) ;
71      }              Paso_SystemMatrix_startRowCollect(A,mask_row) ;
72        }              Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->mainBlock,mask_row,mask_col,main_diagonal_value);
73      }              remote_values=Paso_SystemMatrix_finishColCollect(A);
74    } else {              Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->col_coupleBlock,mask_row,remote_values,0.);
75      if (A->type & MATRIX_FORMAT_CSC) {              remote_values=Paso_SystemMatrix_finishRowCollect(A);
76        #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)              Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->row_coupleBlock,remote_values,mask_col,0.);
77        for (ic=0;ic< A->pattern->n_ptr;ic++) {           }
78      for (iptr=A->pattern->ptr[ic]-index_offset;iptr<A->pattern->ptr[ic+1]-index_offset; iptr++) {         }
79        for (irb=0;irb< A->row_block_size;irb++) {       }
80          irow=irb+A->row_block_size*(A->pattern->index[iptr]-index_offset);    } else {
81          for (icb=0;icb< A->col_block_size;icb++) {       if (A ->col_block_size==1 && A ->row_block_size ==1) {
82            icol=icb+A->col_block_size*ic;         if (A->type & MATRIX_FORMAT_CSC) {
83            if (mask_col[icol]>0. || mask_row[irow]>0. ) {             Paso_SparseMatrix_nullifyRowsAndCols_CSC_BLK1(A->mainBlock,mask_row,mask_col,main_diagonal_value);
84                  l=iptr*A->block_size+irb+A->row_block_size*icb;         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
85          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported.");
86            A->val[l]=main_diagonal_value;         } else {
87          } else {             Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->mainBlock,mask_row,mask_col,main_diagonal_value);
88            A->val[l]=0;         }
89          }       } else {
90            }         if (A->type & MATRIX_FORMAT_CSC) {
91          }             Paso_SparseMatrix_nullifyRowsAndCols_CSC(A->mainBlock,mask_row,mask_col,main_diagonal_value);
92        }         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
93      }             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported with MPI.");
94        }         } else {
95      } else {             Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->mainBlock,mask_row,mask_col,main_diagonal_value);
96        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)         }
97        for (ir=0;ir< A->pattern->n_ptr;ir++) {       }
     for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {  
       for (irb=0;irb< A->row_block_size;irb++) {  
         irow=irb+A->row_block_size*ir;  
         for (icb=0;icb< A->col_block_size;icb++) {  
           icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);  
           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;  
         }  
           }  
         }  
       }  
     }  
       }  
    }  
98    }    }
99    return;    return;
100  }  }
   
 /*  
  * $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.971  
changed lines
  Added in v.1563

  ViewVC Help
Powered by ViewVC 1.1.26