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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 4774 byte(s)
Merging dudley and scons updates from branches

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 Esys_setError(TYPE_ERROR, "Paso_SparseMatrix_invMain: square block size expected.");
46 }
47 if (Esys_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 Esys_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 Esys_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 Esys_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 Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");
121 }
122 */
123 Esys_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

  ViewVC Help
Powered by ViewVC 1.1.26