/[escript]/branches/doubleplusgood/paso/src/SparseMatrix_invMain.cpp
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3120 - (show annotations)
Mon Aug 30 10:48:00 2010 UTC (8 years, 6 months ago) by gross
Original Path: trunk/paso/src/SparseMatrix_invMain.c
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

  ViewVC Help
Powered by ViewVC 1.1.26