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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3482 - (hide 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 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 dhawcroft 631
14 ksteube 1811
15 jgs 150 /**************************************************************/
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 jfenwick 2608 /* Author: Lutz Gross, l.gross@uq.edu.au */
28 jgs 150
29     /**************************************************************/
30    
31     #include "Paso.h"
32     #include "SystemMatrix.h"
33    
34 gross 3369 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 ksteube 1312 register double fac;
82 gross 3369 if (!A->is_balanced) {
83 gross 3303 if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1)) {
84 gross 3369 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 ksteube 1312
94 gross 3369 /* calculat absolute max value over each row */
95 ksteube 1312 #pragma omp parallel for private(irow) schedule(static)
96     for (irow=0; irow<nrow ; ++irow) {
97 gross 3369 A->balance_vector[irow]=0;
98 ksteube 1312 }
99 gross 3369 Paso_SparseMatrix_maxAbsRow_CSR_OFFSET0(A->mainBlock,A->balance_vector);
100 artak 2159 if(A->col_coupleBlock->pattern->ptr!=NULL) {
101 gross 3369 Paso_SparseMatrix_maxAbsRow_CSR_OFFSET0(A->col_coupleBlock,A->balance_vector);
102 artak 2159 }
103 gross 3369
104     /* set balancing vector */
105 ksteube 1312 #pragma omp parallel
106     {
107     #pragma omp for private(irow,fac) schedule(static)
108     for (irow=0; irow<nrow ; ++irow) {
109 gross 3369 fac=A->balance_vector[irow];
110     if (fac>0) {
111     A->balance_vector[irow]=sqrt(1./fac);
112 ksteube 1312 } else {
113 gross 3369 A->balance_vector[irow]=1.;
114 ksteube 1312 }
115     }
116     }
117 gross 3482 {
118 gross 3369 /* 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 jgs 150 }
136 gross 415 }
137 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26