paso/src/SystemMatrix_balancing.c


/*******************************************************
*
* Copyright (c) 2003-2010 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


/**************************************************************/

/* Paso: SystemMatrix                                           */

/*  returns a borrowed reference to the matrix normaliztion vector */

/*  if the vector is in valid a new vector is calculate as the inverse */
/*  of the sum of the absolute value in each row/column                */

/**************************************************************/

/* Copyrights by ACcESS Australia 2005 */
/* Author: Lutz Gross, l.gross@uq.edu.au */

/**************************************************************/

#include "Paso.h"
#include "SystemMatrix.h"

void Paso_SystemMatrix_applyBalanceInPlace(const Paso_SystemMatrix* A, double* x, const bool_t RHS)
{
   const dim_t nrow=A->mainBlock->numRows*A->row_block_size;
   const dim_t ncol=A->mainBlock->numCols*A->col_block_size;
   index_t i;
   
   if (A->is_balanced) {
      if (RHS) {
     #pragma omp parallel for private(i) schedule(static)
     for (i=0; i<nrow ; ++i) {
        x[i]*=A->balance_vector[i];
     }
      } else {
     #pragma omp parallel for private(i) schedule(static)
     for (i=0; i<ncol ; ++i) {
        x[i]*=A->balance_vector[i];
     }
      }
   }
   return;
}

void Paso_SystemMatrix_applyBalance(const Paso_SystemMatrix* A, double* x_out, const double* x, const bool_t RHS)
{
   const dim_t nrow=A->mainBlock->numRows*A->row_block_size;
   const dim_t ncol=A->mainBlock->numCols*A->col_block_size;
   index_t i;
   
   if (A->is_balanced) {
      if (RHS) {
     #pragma omp parallel for private(i) schedule(static)
     for (i=0; i<nrow ; ++i) {
        x_out[i]=x[i]*A->balance_vector[i];
     }
      } else {
     #pragma omp parallel for private(i) schedule(static)
     for (i=0; i<ncol ; ++i) {
        x_out[i]=x[i]*A->balance_vector[i];
     }
      }
   }
   return;
}

void Paso_SystemMatrix_balance(Paso_SystemMatrix* A) {
   dim_t irow;
   const dim_t nrow=A->mainBlock->numRows*A->row_block_size;
   register double fac;
   if (!A->is_balanced) {
      if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1)) {
        Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_balance: No normalization available for compressed sparse column or index offset 1.");
      } 
      if (Esys_checkPtr(A->balance_vector)) {
     Esys_setError(SYSTEM_ERROR,"Paso_SystemMatrix_balance: no memory alloced for balance vector.");
      } 
      if ( ! ( (Paso_SystemMatrix_getGlobalNumRows(A) == Paso_SystemMatrix_getGlobalNumCols(A)) && (A->row_block_size == A->col_block_size) ) ) {
     Esys_setError(SYSTEM_ERROR,"Paso_SystemMatrix_balance: matrix needs to ba a square matrix.");
      }
      if (Esys_noError() ) {

             /* calculat absolute max value over each row */
             #pragma omp parallel for private(irow) schedule(static)
             for (irow=0; irow<nrow ; ++irow) {
        A->balance_vector[irow]=0;
             }
             Paso_SparseMatrix_maxAbsRow_CSR_OFFSET0(A->mainBlock,A->balance_vector);
             if(A->col_coupleBlock->pattern->ptr!=NULL) {
        Paso_SparseMatrix_maxAbsRow_CSR_OFFSET0(A->col_coupleBlock,A->balance_vector);  
             }
             
             /* set balancing vector */
             #pragma omp parallel
             {
                #pragma omp for private(irow,fac) schedule(static)
                for (irow=0; irow<nrow ; ++irow) {
            fac=A->balance_vector[irow];
                    if (fac>0) {
               A->balance_vector[irow]=sqrt(1./fac);
                    } else {
               A->balance_vector[irow]=1.;
                    }
                }
             }
             {
          /* rescale matrix: */
          double *remote_values=NULL;
          /* start exchange */
          Paso_SystemMatrix_startCollect(A, A->balance_vector);
          /* process main block */
          Paso_SparseMatrix_applyDiagonal_CSR_OFFSET0(A->mainBlock,A->balance_vector, A->balance_vector);
          /* finish exchange */
          remote_values=Paso_SystemMatrix_finishCollect(A);
          /* process couple block */
          if (A->col_coupleBlock->pattern->ptr!=NULL) {
             Paso_SparseMatrix_applyDiagonal_CSR_OFFSET0(A->col_coupleBlock,A->balance_vector, remote_values); 
          }
          if (A->row_coupleBlock->pattern->ptr!=NULL) {
             Paso_SparseMatrix_applyDiagonal_CSR_OFFSET0(A->row_coupleBlock, remote_values, A->balance_vector); 
          }
         }
             A->is_balanced=TRUE;
          }
      }
}
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	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision