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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3120 - (hide annotations)
Mon Aug 30 10:48:00 2010 UTC (9 years, 1 month ago) by gross
File MIME type: text/plain
File size: 4727 byte(s)
first iteration on Paso code clean up
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 3094
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 gross 3120 Paso_BlockOps_allMV(n_block,n,block_diag,pivot,x,b);
134 gross 3094 }
135    

  ViewVC Help
Powered by ViewVC 1.1.26