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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3158 - (hide annotations)
Mon Sep 6 06:09:11 2010 UTC (8 years, 7 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 gross 3094
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 gross 3120 #include "BlockOps.h"
30 gross 3158 #include "PasoUtil.h"
31 gross 3094
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 gross 3158 Paso_Copy(n_block*n, x,b);
134     Paso_BlockOps_allMV(n_block,n,block_diag,pivot,x);
135 gross 3094 }
136    

  ViewVC Help
Powered by ViewVC 1.1.26