# Contents of /trunk/paso/src/SparseMatrix_invMain.c

Revision 3094 - (show annotations)
Fri Aug 13 08:38:06 2010 UTC (9 years ago) by gross
File MIME type: text/plain
File size: 4722 byte(s)
```The MPI and sequational GAUSS_SEIDEL have been merged.
The couring and main diagonal pointer is now manged by the patternm which means that they are calculated once only even if the preconditioner is deleted.

```
 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: inverts the main diagonal of a SparseMatrix: */ 18 19 /**************************************************************/ 20 21 /* Copyrights by ACcESS Australia 2010 */ 22 /* Author: Lutz Gross, l.gross@uq.edu.au */ 23 24 /**************************************************************/ 25 26 #include "Paso.h" 27 #include "SparseMatrix.h" 28 #include "Solver.h" 29 30 void Paso_SparseMatrix_invMain(Paso_SparseMatrix * A_p, double* inv_diag, int* pivot) { 31 /* inv_diag=MEMALLOC( A->numRows * A_p-> block_size,double); 32 pivot=MEMALLOC( A->numRows * A->row_block_size */ 33 dim_t n=A_p->numRows; 34 dim_t n_block=A_p->row_block_size; 35 dim_t m_block=A_p->col_block_size; 36 /* dim_t block_size=A_p->block_size; */ 37 dim_t i; 38 register index_t iPtr; 39 register double A11,A12,A13,A21,A22,A23,A31,A32,A33,D; 40 index_t* main_ptr=Paso_Pattern_borrowMainDiagonalPointer(A_p->pattern); 41 /* check matrix is square */ 42 if (m_block != n_block) { 43 Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_invMain: square block size expected."); 44 } 45 if (Paso_noError()) { 46 47 if (n_block==1) { 48 #pragma omp parallel for private(i, iPtr, A11) schedule(static) 49 for (i = 0; i < n; i++) { 50 iPtr= main_ptr[i]; 51 A11=A_p->val[iPtr]; 52 if ( ABS(A11) > 0.) { 53 inv_diag[i]=1./A11; 54 } else { 55 Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block."); 56 } 57 } 58 } else if (n_block==2) { 59 #pragma omp parallel for private(i, iPtr, A11,A12,A21,A22,D) schedule(static) 60 for (i = 0; i < n; i++) { 61 iPtr= main_ptr[i]; 62 A11=A_p->val[iPtr*4]; 63 A12=A_p->val[iPtr*4+2]; 64 A21=A_p->val[iPtr*4+1]; 65 A22=A_p->val[iPtr*4+3]; 66 D = A11*A22-A12*A21; 67 if (ABS(D) > 0 ){ 68 D=1./D; 69 inv_diag[i*4]= A22*D; 70 inv_diag[i*4+1]=-A21*D; 71 inv_diag[i*4+2]=-A12*D; 72 inv_diag[i*4+3]= A11*D; 73 } else { 74 Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block."); 75 } 76 } 77 } else if (n_block==3) { 78 #pragma omp parallel for private(i, iPtr,A11,A12,A13,A21,A22,A23,A31,A32,A33,D) schedule(static) 79 for (i = 0; i < n; i++) { 80 iPtr= main_ptr[i]; 81 A11=A_p->val[iPtr*9 ]; 82 A21=A_p->val[iPtr*9+1]; 83 A31=A_p->val[iPtr*9+2]; 84 A12=A_p->val[iPtr*9+3]; 85 A22=A_p->val[iPtr*9+4]; 86 A32=A_p->val[iPtr*9+5]; 87 A13=A_p->val[iPtr*9+6]; 88 A23=A_p->val[iPtr*9+7]; 89 A33=A_p->val[iPtr*9+8]; 90 D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22); 91 if (ABS(D) > 0 ){ 92 D=1./D; 93 inv_diag[i*9 ]=(A22*A33-A23*A32)*D; 94 inv_diag[i*9+1]=(A31*A23-A21*A33)*D; 95 inv_diag[i*9+2]=(A21*A32-A31*A22)*D; 96 inv_diag[i*9+3]=(A13*A32-A12*A33)*D; 97 inv_diag[i*9+4]=(A11*A33-A31*A13)*D; 98 inv_diag[i*9+5]=(A12*A31-A11*A32)*D; 99 inv_diag[i*9+6]=(A12*A23-A13*A22)*D; 100 inv_diag[i*9+7]=(A13*A21-A11*A23)*D; 101 inv_diag[i*9+8]=(A11*A22-A12*A21)*D; 102 } else { 103 Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block."); 104 } 105 } 106 } else { 107 /* 108 #pragma omp parallel for private(i, INFO) schedule(static) 109 for (i = 0; i < n; i++) { 110 iPtr= main_ptr[i]; 111 memcpy((void*)&(diag[i*block_size], (void*)(&A_p->val[iPtr*block_size]), sizeof(double)*(size_t)block_size); 112 113 int dgetrf_(integer *m, integer *n, doublereal *a, integer * 114 lda, integer *ipiv, integer *info); 115 116 dgetrf_(n_block, m_block, &(diag[i*block_size], n_block, pivot[i*n_block], info ); 117 if (INFO >0 ) { 118 Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block."); 119 } 120 */ 121 Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_invMain: Right now there is support block size less than 4 only"); 122 } 123 } 124 } 125 void Paso_SparseMatrix_applyBlockMatrix(Paso_SparseMatrix * A_p, double* block_diag, int* pivot, double*x, double *b) { 126 127 /* inv_diag=MEMALLOC( A->numRows * A_p-> block_size,double); 128 pivot=MEMALLOC( A->numRows * A->row_block_size */ 129 dim_t n=A_p->numRows; 130 dim_t n_block=A_p->row_block_size; 131 132 Paso_Solver_applyBlockDiagonalMatrix(n_block,n,block_diag,pivot,x,b); 133 } 134