/[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 3094 - (hide annotations)
Fri Aug 13 08:38:06 2010 UTC (9 years ago) by gross
File MIME type: text/plain
File size: 4722 byte(s)
The MPI and sequational GAUSS_SEIDEL have been merged.
The couring and main diagonal pointer is now manged by the patternm which means that they are calculated once only even if the preconditioner is deleted.



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

  ViewVC Help
Powered by ViewVC 1.1.26