Contents of /trunk/paso/src/SystemMatrix_balancing.c

Revision 3482 - (show annotations)
Wed Mar 23 04:06:52 2011 UTC (8 years, 6 months ago) by gross
File MIME type: text/plain
File size: 4732 byte(s)
```some work on AMG MPI
```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2010 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 */ 18 19 /* returns a borrowed reference to the matrix normaliztion vector */ 20 21 /* if the vector is in valid a new vector is calculate as the inverse */ 22 /* of the sum of the absolute value in each row/column */ 23 24 /**************************************************************/ 25 26 /* Copyrights by ACcESS Australia 2005 */ 27 /* Author: Lutz Gross, l.gross@uq.edu.au */ 28 29 /**************************************************************/ 30 31 #include "Paso.h" 32 #include "SystemMatrix.h" 33 34 void Paso_SystemMatrix_applyBalanceInPlace(const Paso_SystemMatrix* A, double* x, const bool_t RHS) 35 { 36 const dim_t nrow=A->mainBlock->numRows*A->row_block_size; 37 const dim_t ncol=A->mainBlock->numCols*A->col_block_size; 38 index_t i; 39 40 if (A->is_balanced) { 41 if (RHS) { 42 #pragma omp parallel for private(i) schedule(static) 43 for (i=0; ibalance_vector[i]; 45 } 46 } else { 47 #pragma omp parallel for private(i) schedule(static) 48 for (i=0; ibalance_vector[i]; 50 } 51 } 52 } 53 return; 54 } 55 56 void Paso_SystemMatrix_applyBalance(const Paso_SystemMatrix* A, double* x_out, const double* x, const bool_t RHS) 57 { 58 const dim_t nrow=A->mainBlock->numRows*A->row_block_size; 59 const dim_t ncol=A->mainBlock->numCols*A->col_block_size; 60 index_t i; 61 62 if (A->is_balanced) { 63 if (RHS) { 64 #pragma omp parallel for private(i) schedule(static) 65 for (i=0; ibalance_vector[i]; 67 } 68 } else { 69 #pragma omp parallel for private(i) schedule(static) 70 for (i=0; ibalance_vector[i]; 72 } 73 } 74 } 75 return; 76 } 77 78 void Paso_SystemMatrix_balance(Paso_SystemMatrix* A) { 79 dim_t irow; 80 const dim_t nrow=A->mainBlock->numRows*A->row_block_size; 81 register double fac; 82 if (!A->is_balanced) { 83 if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1)) { 84 Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_balance: No normalization available for compressed sparse column or index offset 1."); 85 } 86 if (Esys_checkPtr(A->balance_vector)) { 87 Esys_setError(SYSTEM_ERROR,"Paso_SystemMatrix_balance: no memory alloced for balance vector."); 88 } 89 if ( ! ( (Paso_SystemMatrix_getGlobalNumRows(A) == Paso_SystemMatrix_getGlobalNumCols(A)) && (A->row_block_size == A->col_block_size) ) ) { 90 Esys_setError(SYSTEM_ERROR,"Paso_SystemMatrix_balance: matrix needs to ba a square matrix."); 91 } 92 if (Esys_noError() ) { 93 94 /* calculat absolute max value over each row */ 95 #pragma omp parallel for private(irow) schedule(static) 96 for (irow=0; irowbalance_vector[irow]=0; 98 } 99 Paso_SparseMatrix_maxAbsRow_CSR_OFFSET0(A->mainBlock,A->balance_vector); 100 if(A->col_coupleBlock->pattern->ptr!=NULL) { 101 Paso_SparseMatrix_maxAbsRow_CSR_OFFSET0(A->col_coupleBlock,A->balance_vector); 102 } 103 104 /* set balancing vector */ 105 #pragma omp parallel 106 { 107 #pragma omp for private(irow,fac) schedule(static) 108 for (irow=0; irowbalance_vector[irow]; 110 if (fac>0) { 111 A->balance_vector[irow]=sqrt(1./fac); 112 } else { 113 A->balance_vector[irow]=1.; 114 } 115 } 116 } 117 { 118 /* rescale matrix: */ 119 double *remote_values=NULL; 120 /* start exchange */ 121 Paso_SystemMatrix_startCollect(A, A->balance_vector); 122 /* process main block */ 123 Paso_SparseMatrix_applyDiagonal_CSR_OFFSET0(A->mainBlock,A->balance_vector, A->balance_vector); 124 /* finish exchange */ 125 remote_values=Paso_SystemMatrix_finishCollect(A); 126 /* process couple block */ 127 if (A->col_coupleBlock->pattern->ptr!=NULL) { 128 Paso_SparseMatrix_applyDiagonal_CSR_OFFSET0(A->col_coupleBlock,A->balance_vector, remote_values); 129 } 130 if (A->row_coupleBlock->pattern->ptr!=NULL) { 131 Paso_SparseMatrix_applyDiagonal_CSR_OFFSET0(A->row_coupleBlock, remote_values, A->balance_vector); 132 } 133 } 134 A->is_balanced=TRUE; 135 } 136 } 137 }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision