/[escript]/trunk/paso/src/SystemMatrix_balancing.cpp
ViewVC logotype

Diff of /trunk/paso/src/SystemMatrix_balancing.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4828 by caltinay, Thu Mar 27 07:00:30 2014 UTC revision 4829 by caltinay, Thu Apr 3 04:02:53 2014 UTC
# Line 15  Line 15 
15  *****************************************************************************/  *****************************************************************************/
16    
17    
18  /************************************************************************************/  /****************************************************************************/
19    
20  /* Paso: SystemMatrix                                         */  /* Paso: SystemMatrix                                         */
21    
# Line 24  Line 24 
24  /*  If the vector is in valid a new vector is calculated as the  */  /*  If the vector is in valid a new vector is calculated as the  */
25  /*  inverse of the sum of the absolute value in each row/column. */  /*  inverse of the sum of the absolute value in each row/column. */
26    
27  /************************************************************************************/  /****************************************************************************/
28    
29  /* Copyrights by ACcESS Australia 2005 */  /* Copyrights by ACcESS Australia 2005 */
30  /* Author: Lutz Gross, l.gross@uq.edu.au */  /* Author: Lutz Gross, l.gross@uq.edu.au */
31    
32  /************************************************************************************/  /****************************************************************************/
33    
34  #include "Paso.h"  #include "Paso.h"
35  #include "SystemMatrix.h"  #include "SystemMatrix.h"
# Line 42  void Paso_SystemMatrix_applyBalanceInPla Line 42  void Paso_SystemMatrix_applyBalanceInPla
42        
43     if (A->is_balanced) {     if (A->is_balanced) {
44        if (RHS) {        if (RHS) {
45       #pragma omp parallel for private(i) schedule(static)           #pragma omp parallel for private(i) schedule(static)
46       for (i=0; i<nrow ; ++i) {           for (i=0; i<nrow ; ++i) {
47          x[i]*=A->balance_vector[i];              x[i]*=A->balance_vector[i];
48       }           }
49        } else {        } else {
50       #pragma omp parallel for private(i) schedule(static)           #pragma omp parallel for private(i) schedule(static)
51       for (i=0; i<ncol ; ++i) {           for (i=0; i<ncol ; ++i) {
52          x[i]*=A->balance_vector[i];              x[i]*=A->balance_vector[i];
53       }           }
54        }        }
55     }     }
56     return;     return;
# Line 64  void Paso_SystemMatrix_applyBalance(cons Line 64  void Paso_SystemMatrix_applyBalance(cons
64        
65     if (A->is_balanced) {     if (A->is_balanced) {
66        if (RHS) {        if (RHS) {
67       #pragma omp parallel for private(i) schedule(static)           #pragma omp parallel for private(i) schedule(static)
68       for (i=0; i<nrow ; ++i) {           for (i=0; i<nrow ; ++i) {
69          x_out[i]=x[i]*A->balance_vector[i];              x_out[i]=x[i]*A->balance_vector[i];
70       }           }
71        } else {        } else {
72       #pragma omp parallel for private(i) schedule(static)           #pragma omp parallel for private(i) schedule(static)
73       for (i=0; i<ncol ; ++i) {           for (i=0; i<ncol ; ++i) {
74          x_out[i]=x[i]*A->balance_vector[i];              x_out[i]=x[i]*A->balance_vector[i];
75       }           }
76        }        }
77     }     }
78     return;     return;
79  }  }
80    
81  void Paso_SystemMatrix_balance(Paso_SystemMatrix* A) {  void Paso_SystemMatrix_balance(Paso_SystemMatrix* A)
82     dim_t irow;  {
83     const dim_t nrow=Paso_SystemMatrix_getTotalNumRows(A);      dim_t irow;
84     register double fac;      const dim_t nrow=Paso_SystemMatrix_getTotalNumRows(A);
85        register double fac;
86     if (!A->is_balanced) {  
87        if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1)) {      if (!A->is_balanced) {
88          Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_balance: No normalization available for compressed sparse column or index offset 1.");          if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1)) {
89        }              Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_balance: No normalization available for compressed sparse column or index offset 1.");
90        if ( ! ( (Paso_SystemMatrix_getGlobalNumRows(A) == Paso_SystemMatrix_getGlobalNumCols(A)) && (A->row_block_size == A->col_block_size) ) ) {          }
91       Esys_setError(SYSTEM_ERROR,"Paso_SystemMatrix_balance: matrix needs to be a square matrix.");          if ( ! ( (Paso_SystemMatrix_getGlobalNumRows(A) == Paso_SystemMatrix_getGlobalNumCols(A)) && (A->row_block_size == A->col_block_size) ) ) {
92        }              Esys_setError(SYSTEM_ERROR,"Paso_SystemMatrix_balance: matrix needs to be a square matrix.");
93        if (Esys_noError() ) {          }
94            if (Esys_noError() ) {
95               /* calculate absolute max value over each row */              /* calculate absolute max value over each row */
96               #pragma omp parallel for private(irow) schedule(static)              #pragma omp parallel for private(irow) schedule(static)
97               for (irow=0; irow<nrow ; ++irow) {              for (irow=0; irow<nrow ; ++irow) {
98          A->balance_vector[irow]=0;                  A->balance_vector[irow]=0;
99               }              }
100               paso::SparseMatrix_maxAbsRow_CSR_OFFSET0(A->mainBlock,A->balance_vector);              A->mainBlock->maxAbsRow_CSR_OFFSET0(A->balance_vector);
101               if(A->col_coupleBlock->pattern->ptr!=NULL) {              if(A->col_coupleBlock->pattern->ptr!=NULL) {
102                   paso::SparseMatrix_maxAbsRow_CSR_OFFSET0(A->col_coupleBlock,A->balance_vector);                     A->col_coupleBlock->maxAbsRow_CSR_OFFSET0(A->balance_vector);
103               }              }
104                            
105               /* set balancing vector */              /* set balancing vector */
106               #pragma omp parallel              #pragma omp parallel
107               {              {
108                  #pragma omp for private(irow,fac) schedule(static)                  #pragma omp for private(irow,fac) schedule(static)
109                  for (irow=0; irow<nrow ; ++irow) {                  for (irow=0; irow<nrow ; ++irow) {
110              fac=A->balance_vector[irow];                      fac=A->balance_vector[irow];
111                      if (fac>0) {                      if (fac>0) {
112                 A->balance_vector[irow]=sqrt(1./fac);                         A->balance_vector[irow]=sqrt(1./fac);
113                      } else {                      } else {
114                 A->balance_vector[irow]=1.;                         A->balance_vector[irow]=1.;
115                      }                      }
116                  }                  }
117          }              }
118              {              {
119            /* rescale matrix: */                    /* rescale matrix: */
120            double *remote_values=NULL;                    double *remote_values=NULL;
121            /* start exchange */                    /* start exchange */
122            Paso_SystemMatrix_startCollect(A, A->balance_vector);                    Paso_SystemMatrix_startCollect(A, A->balance_vector);
123            /* process main block */                    /* process main block */
124            paso::SparseMatrix_applyDiagonal_CSR_OFFSET0(A->mainBlock,A->balance_vector, A->balance_vector);                    A->mainBlock->applyDiagonal_CSR_OFFSET0(A->balance_vector, A->balance_vector);
125            /* finish exchange */                    /* finish exchange */
126            remote_values=Paso_SystemMatrix_finishCollect(A);                    remote_values=Paso_SystemMatrix_finishCollect(A);
127            /* process couple block */                    /* process couple block */
128            if (A->col_coupleBlock->pattern->ptr!=NULL) {                    if (A->col_coupleBlock->pattern->ptr!=NULL) {
129                paso::SparseMatrix_applyDiagonal_CSR_OFFSET0(A->col_coupleBlock,A->balance_vector, remote_values);                      A->col_coupleBlock->applyDiagonal_CSR_OFFSET0(A->balance_vector, remote_values);
130            }                    }
131            if (A->row_coupleBlock->pattern->ptr!=NULL) {                    if (A->row_coupleBlock->pattern->ptr!=NULL) {
132                paso::SparseMatrix_applyDiagonal_CSR_OFFSET0(A->row_coupleBlock, remote_values, A->balance_vector);                      A->row_coupleBlock->applyDiagonal_CSR_OFFSET0(remote_values, A->balance_vector);
133            }                    }
134           }              }
135               A->is_balanced=TRUE;              A->is_balanced=true;
136            }          }
137        }      }
138  }  }
139    

Legend:
Removed from v.4828  
changed lines
  Added in v.4829

  ViewVC Help
Powered by ViewVC 1.1.26