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

Diff of /trunk/paso/src/UMFPACK.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3009 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 3010 by artak, Tue Apr 27 05:10:46 2010 UTC
# Line 30  Line 30 
30  #include <omp.h>  #include <omp.h>
31  #endif  #endif
32    
33    #define MAX_BLOCK_SIZE 3
34    
35  /**************************************************************/  /**************************************************************/
36    
37  /*  free any extra stuff possibly used by the UMFPACK library */  /*  free any extra stuff possibly used by the UMFPACK library */
# Line 93  void Paso_UMFPACK(Paso_SystemMatrix* A, Line 95  void Paso_UMFPACK(Paso_SystemMatrix* A,
95  void Paso_UMFPACK1(Paso_UMFPACK_Handler** pt, Paso_SparseMatrix* A, double* out, double* in, const int refines) {  void Paso_UMFPACK1(Paso_UMFPACK_Handler** pt, Paso_SparseMatrix* A, double* out, double* in, const int refines) {
96       double control[UMFPACK_CONTROL], info[UMFPACK_INFO];       double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
97       int error = UMFPACK_OK;       int error = UMFPACK_OK;
98       umfpack_di_defaults(control);       dim_t i,j;
99         /*****/
100         double **xx;
101         double **bb;
102         Paso_SparseMatrix *block[MAX_BLOCK_SIZE];
103        
104         xx=MEMALLOC(A->row_block_size,double*);
105         bb=MEMALLOC(A->row_block_size,double*);
106         if (Paso_checkPtr(xx) || Paso_checkPtr(bb)) return;
107         /****/
108    
109       if (*pt==NULL) {       umfpack_di_defaults(control);
110          int n = A->numRows;      
111          *pt=(MEMALLOC(1,Paso_UMFPACK_Handler));       for (i=0;i<A->row_block_size;i++) {
112          if (Paso_checkPtr(*pt)) return;            xx[i]=MEMALLOC(A->numRows,double);
113          /* call LDU symbolic factorization: */            bb[i]=MEMALLOC(A->numRows,double);
114          error=umfpack_di_symbolic(n,n,A->pattern->ptr,A->pattern->index,A->val,&((*pt)->symbolic),control,info);            if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
115          if (error != UMFPACK_OK) {       }
116               Paso_setError(VALUE_ERROR,"symbolic factorization failed.");      
117               return;       /*#pragma omp parallel for private(i,j) schedule(static)*/
118          } else {       for (i=0;i<A->numRows;i++) {
119              /* call LDU factorization: */           for (j=0;j<A->row_block_size;j++) {
120              error= umfpack_di_numeric(A->pattern->ptr,A->pattern->index,A->val,(*pt)->symbolic,&((*pt)->numeric),control,info);            bb[j][i]=in[A->row_block_size*i+j];
121             if (error != UMFPACK_OK) {            xx[j][i]=0;  
122               Paso_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");           }
123               return;        }
124             }      
125          }       for (i=0;i<MAX_BLOCK_SIZE;++i) {
126       }            block[i]=NULL;
127       if (Paso_noError())  {       }
128          /* call forward backward substitution: */      
129          control[UMFPACK_IRSTEP]=refines; /* number of refinement steps */       for (i=0;i<A->row_block_size;++i) {
130          error=umfpack_di_solve(UMFPACK_A,A->pattern->ptr,A->pattern->index,A->val,out,in,(*pt)->numeric,control,info);          block[i]=Paso_SparseMatrix_getBlock(A,i+1);
131          if (error != UMFPACK_OK) {  
132                Paso_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");            if (*pt==NULL) {
133                return;               int n = block[i]->numRows;
134          }               *pt=(MEMALLOC(1,Paso_UMFPACK_Handler));
135                 if (Paso_checkPtr(*pt)) return;
136                 /* call LDU symbolic factorization: */
137                 error=umfpack_di_symbolic(n,n,block[i]->pattern->ptr,block[i]->pattern->index,block[i]->val,&((*pt)->symbolic),control,info);
138                 if (error != UMFPACK_OK) {
139                      Paso_setError(VALUE_ERROR,"symbolic factorization failed.");
140                      return;
141                 } else {
142                     /* call LDU factorization: */
143                     error= umfpack_di_numeric(block[i]->pattern->ptr,block[i]->pattern->index,block[i]->val,(*pt)->symbolic,&((*pt)->numeric),control,info);
144                    if (error != UMFPACK_OK) {
145                      Paso_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");
146                      return;
147                    }
148                 }
149              }
150              
151              if (Paso_noError())  {
152                 /* call forward backward substitution: */
153                 control[UMFPACK_IRSTEP]=refines; /* number of refinement steps */
154                 error=umfpack_di_solve(UMFPACK_A,block[i]->pattern->ptr,block[i]->pattern->index,block[i]->val,xx[i],bb[i],(*pt)->numeric,control,info);
155                 if (error != UMFPACK_OK) {
156                       Paso_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");
157                       return;
158                 }
159              }
160        
161       }       }
162        
163         /*#pragma omp parallel for private(i,j) schedule(static)*/
164          for (i=0;i<A->numRows;i++) {
165              for (j=0;j<A->row_block_size;j++) {
166              out[A->row_block_size*i+j]=xx[j][i];
167              }
168           }
169          
170         for (i=0;i<A->row_block_size;i++) {
171                    MEMFREE(xx[i]);
172                    MEMFREE(bb[i]);
173                    Paso_SparseMatrix_free(block[i]);
174          }
175          
176  #else  #else
177      Paso_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK is not avialble.");      Paso_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK is not avialble.");
178  #endif  #endif
Legend:
Removed from v.3009  
changed lines
  Added in v.3010
  ViewVC Help
Powered by ViewVC 1.1.26