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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3369 - (show annotations)
Fri Nov 19 06:26:11 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/plain
File size: 4733 byte(s)
matrix balance introduced in paso
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; i<nrow ; ++i) {
44 x[i]*=A->balance_vector[i];
45 }
46 } else {
47 #pragma omp parallel for private(i) schedule(static)
48 for (i=0; i<ncol ; ++i) {
49 x[i]*=A->balance_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; i<nrow ; ++i) {
66 x_out[i]=x[i]*A->balance_vector[i];
67 }
68 } else {
69 #pragma omp parallel for private(i) schedule(static)
70 for (i=0; i<ncol ; ++i) {
71 x_out[i]=x[i]*A->balance_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; irow<nrow ; ++irow) {
97 A->balance_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; irow<nrow ; ++irow) {
109 fac=A->balance_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

  ViewVC Help
Powered by ViewVC 1.1.26