# Contents of /branches/doubleplusgood/paso/src/SparseMatrix_invMain.cpp

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