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

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