/[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 temp_trunk_copy/paso/src/SystemMatrix_nullifyRowsAndCols.c revision 1384 by phornby, Fri Jan 11 02:29:38 2008 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /* $Id$ */
3    
4    /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  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                                       */
# Line 24  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 ir,icol,iptr,icb,irb,irow,ic,l;    double *remote_values=NULL;
41      Paso_MPIInfo *mpi_info=A->mpi_info;
42    if (A ->col_block_size==1 && A ->row_block_size ==1) {    if (mpi_info->size>1) {
43      switch(A->type) {       if (A ->col_block_size==1 && A ->row_block_size ==1) {
44      case CSR:         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 (irow=0;irow< A->pattern->n_ptr;irow++) {             return;
47      /* TODO: test mask_row here amd not inside every loop */         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
48      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.");
49        icol=A->pattern->index[iptr]-INDEX_OFFSET;             return;
50        if (mask_col[icol]>0. || mask_row[irow]>0. ) {         } else {
51              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {           Paso_SystemMatrix_allocBuffer(A);
52            A->val[iptr]=main_diagonal_value;           if (Paso_noError()) {
53              } else {              Paso_SystemMatrix_startCollect(A,mask_col) ;
54            A->val[iptr]=0;              Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->mainBlock,mask_row,mask_col,main_diagonal_value);
55              }              remote_values=Paso_SystemMatrix_finishCollect(A);
56        }              Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->coupleBlock,mask_row,remote_values,0.);
57      }           }
58        }           Paso_SystemMatrix_freeBuffer(A);
59        break;         }
60      case CSC:       } else {
61        #pragma omp parallel for private(irow, iptr,icol) schedule(static)         if (A->type & MATRIX_FORMAT_CSC) {
62        for (icol=0;icol< A->pattern->n_ptr;icol++) {             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: CSC is not supported by MPI.");
63      for (iptr=A->pattern->ptr[icol]-PTR_OFFSET;iptr<A->pattern->ptr[icol+1]-PTR_OFFSET; iptr++) {             return;
64        irow=A->pattern->index[iptr]-INDEX_OFFSET;         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
65        if (mask_col[icol]>0. || mask_row[irow]>0. ) {             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported with MPI.");
66              if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {             return;
67            A->val[iptr]=main_diagonal_value;         } else {
68              } else {           Paso_SystemMatrix_allocBuffer(A);
69            A->val[iptr]=0;           if (Paso_noError()) {
70              }              Paso_SystemMatrix_startCollect(A,mask_col) ;
71        }              Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->mainBlock,mask_row,mask_col,main_diagonal_value);
72      }              remote_values=Paso_SystemMatrix_finishCollect(A);
73        }              Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->coupleBlock,mask_row,remote_values,0.);
74        break;           }
75      default:           Paso_SystemMatrix_freeBuffer(A);
76        Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");         }
77      } /* switch A->type */       }
78    } else {    } else {
79      switch(A->type) {       if (A ->col_block_size==1 && A ->row_block_size ==1) {
80      case CSR:         if (A->type & MATRIX_FORMAT_CSC) {
81        #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)             Paso_SparseMatrix_nullifyRowsAndCols_CSC_BLK1(A->mainBlock,mask_row,mask_col,main_diagonal_value);
82        for (ir=0;ir< A->pattern->n_ptr;ir++) {         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
83      for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported.");
84        for (irb=0;irb< A->row_block_size;irb++) {         } else {
85          irow=irb+A->row_block_size*ir;             Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(A->mainBlock,mask_row,mask_col,main_diagonal_value);
86          for (icb=0;icb< A->col_block_size;icb++) {         }
87            icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);       } else {
88            if (mask_col[icol]>0. || mask_row[irow]>0. ) {         if (A->type & MATRIX_FORMAT_CSC) {
89                  l=iptr*A->block_size+irb+A->row_block_size*icb;             Paso_SparseMatrix_nullifyRowsAndCols_CSC(A->mainBlock,mask_row,mask_col,main_diagonal_value);
90          if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {         } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
91            A->val[l]=main_diagonal_value;             Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_nullifyRowsAndCols: TRILINOS is not supported with MPI.");
92          } else {         } else {
93            A->val[l]=0;             Paso_SparseMatrix_nullifyRowsAndCols_CSR(A->mainBlock,mask_row,mask_col,main_diagonal_value);
94          }         }
95            }       }
         }  
       }  
     }  
       }  
       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 */  
96    }    }
97    return;    return;
98  }  }
   
 /*  
  * $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.1384

  ViewVC Help
Powered by ViewVC 1.1.26