/[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 3642 - (show annotations)
Thu Oct 27 03:41:51 2011 UTC (7 years, 5 months ago) by caltinay
File MIME type: text/plain
File size: 4701 byte(s)
Assorted spelling/comment fixes 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 normalization vector. */
20
21 /* If the vector is in valid a new vector is calculated as the */
22 /* inverse 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=Paso_SystemMatrix_getTotalNumRows(A);
37 const dim_t ncol=Paso_SystemMatrix_getTotalNumCols(A);
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=Paso_SystemMatrix_getTotalNumRows(A);
59 const dim_t ncol=Paso_SystemMatrix_getTotalNumCols(A);
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=Paso_SystemMatrix_getTotalNumRows(A);
81 register double fac;
82
83 if (!A->is_balanced) {
84 if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1)) {
85 Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_balance: No normalization available for compressed sparse column or index offset 1.");
86 }
87 if (Esys_checkPtr(A->balance_vector)) {
88 Esys_setError(SYSTEM_ERROR,"Paso_SystemMatrix_balance: no memory allocated for balance vector.");
89 }
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.");
92 }
93 if (Esys_noError() ) {
94
95 /* calculate absolute max value over each row */
96 #pragma omp parallel for private(irow) schedule(static)
97 for (irow=0; irow<nrow ; ++irow) {
98 A->balance_vector[irow]=0;
99 }
100 Paso_SparseMatrix_maxAbsRow_CSR_OFFSET0(A->mainBlock,A->balance_vector);
101 if(A->col_coupleBlock->pattern->ptr!=NULL) {
102 Paso_SparseMatrix_maxAbsRow_CSR_OFFSET0(A->col_coupleBlock,A->balance_vector);
103 }
104
105 /* set balancing vector */
106 #pragma omp parallel
107 {
108 #pragma omp for private(irow,fac) schedule(static)
109 for (irow=0; irow<nrow ; ++irow) {
110 fac=A->balance_vector[irow];
111 if (fac>0) {
112 A->balance_vector[irow]=sqrt(1./fac);
113 } else {
114 A->balance_vector[irow]=1.;
115 }
116 }
117 }
118 {
119 /* rescale matrix: */
120 double *remote_values=NULL;
121 /* start exchange */
122 Paso_SystemMatrix_startCollect(A, A->balance_vector);
123 /* process main block */
124 Paso_SparseMatrix_applyDiagonal_CSR_OFFSET0(A->mainBlock,A->balance_vector, A->balance_vector);
125 /* finish exchange */
126 remote_values=Paso_SystemMatrix_finishCollect(A);
127 /* process couple block */
128 if (A->col_coupleBlock->pattern->ptr!=NULL) {
129 Paso_SparseMatrix_applyDiagonal_CSR_OFFSET0(A->col_coupleBlock,A->balance_vector, remote_values);
130 }
131 if (A->row_coupleBlock->pattern->ptr!=NULL) {
132 Paso_SparseMatrix_applyDiagonal_CSR_OFFSET0(A->row_coupleBlock, remote_values, A->balance_vector);
133 }
134 }
135 A->is_balanced=TRUE;
136 }
137 }
138 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26