/[escript]/branches/domexper/paso/src/RILU.c
ViewVC logotype

Diff of /branches/domexper/paso/src/RILU.c

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

trunk/paso/src/RILU.c revision 3232 by gross, Mon Sep 6 09:17:28 2010 UTC branches/domexper/paso/src/RILU.c revision 3234 by jfenwick, Mon Oct 4 01:46:30 2010 UTC
# Line 99  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 99  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
99    out->x_C=NULL;    out->x_C=NULL;
100    out->b_C=NULL;    out->b_C=NULL;
101    
102    if ( !(Paso_checkPtr(mis_marker) || Paso_checkPtr(out) || Paso_checkPtr(counter) ) ) {    if ( !(Esys_checkPtr(mis_marker) || Esys_checkPtr(out) || Esys_checkPtr(counter) ) ) {
103       /* identify independend set of rows/columns */       /* identify independend set of rows/columns */
104       time0=Paso_timer();       time0=Esys_timer();
105       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
106       for (i=0;i<n;++i) mis_marker[i]=-1;       for (i=0;i<n;++i) mis_marker[i]=-1;
107       Paso_Pattern_mis(A_p->pattern,mis_marker);       Paso_Pattern_mis(A_p->pattern,mis_marker);
108       time2=Paso_timer()-time0;       time2=Esys_timer()-time0;
109       if (Paso_noError()) {       if (Esys_noError()) {
110          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
111          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
112          out->n=n;          out->n=n;
# Line 116  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 116  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
116          out->rows_in_F=MEMALLOC(out->n_F,index_t);          out->rows_in_F=MEMALLOC(out->n_F,index_t);
117          out->inv_A_FF=MEMALLOC(n_block*n_block*out->n_F,double);          out->inv_A_FF=MEMALLOC(n_block*n_block*out->n_F,double);
118          out->A_FF_pivot=NULL; /* later use for block size>3 */          out->A_FF_pivot=NULL; /* later use for block size>3 */
119          if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->inv_A_FF) || Paso_checkPtr(out->rows_in_F) ) ) {          if (! (Esys_checkPtr(out->mask_F) || Esys_checkPtr(out->inv_A_FF) || Esys_checkPtr(out->rows_in_F) ) ) {
120             #pragma omp parallel             #pragma omp parallel
121             {             {
122                /* creates an index for F from mask */                /* creates an index for F from mask */
# Line 142  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 142  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
142                                          sizeof(index_t),                                          sizeof(index_t),
143                                          Paso_comparIndex);                                          Paso_comparIndex);
144                  if (where_p==NULL) {                  if (where_p==NULL) {
145                      Paso_setError(VALUE_ERROR, "Paso_Solver_getRILU: main diagonal element missing.");                      Esys_setError(VALUE_ERROR, "Paso_Solver_getRILU: main diagonal element missing.");
146                  } else {                  } else {
147                      iPtr+=(index_t)(where_p-index);                      iPtr+=(index_t)(where_p-index);
148                      /* get inverse of A_FF block: */                      /* get inverse of A_FF block: */
# Line 150  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 150  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
150                         if (ABS(A_p->val[iPtr])>0.) {                         if (ABS(A_p->val[iPtr])>0.) {
151                              out->inv_A_FF[i]=1./A_p->val[iPtr];                              out->inv_A_FF[i]=1./A_p->val[iPtr];
152                         } else {                         } else {
153                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getRILU: Break-down in RILU decomposition: non-regular main diagonal block.");                              Esys_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getRILU: Break-down in RILU decomposition: non-regular main diagonal block.");
154                         }                         }
155                      } else if (n_block==2) {                      } else if (n_block==2) {
156                         A11=A_p->val[iPtr*4];                         A11=A_p->val[iPtr*4];
# Line 165  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 165  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
165                              out->inv_A_FF[i*4+2]=-A12*D;                              out->inv_A_FF[i*4+2]=-A12*D;
166                              out->inv_A_FF[i*4+3]= A11*D;                              out->inv_A_FF[i*4+3]= A11*D;
167                         } else {                         } else {
168                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getRILU:Break-down in RILU decomposition: non-regular main diagonal block.");                              Esys_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getRILU:Break-down in RILU decomposition: non-regular main diagonal block.");
169                         }                         }
170                      } else if (n_block==3) {                      } else if (n_block==3) {
171                         A11=A_p->val[iPtr*9  ];                         A11=A_p->val[iPtr*9  ];
# Line 190  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 190  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
190                              out->inv_A_FF[i*9+7]=(A13*A21-A11*A23)*D;                              out->inv_A_FF[i*9+7]=(A13*A21-A11*A23)*D;
191                              out->inv_A_FF[i*9+8]=(A11*A22-A12*A21)*D;                              out->inv_A_FF[i*9+8]=(A11*A22-A12*A21)*D;
192                         } else {                         } else {
193                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getRILU:Break-down in RILU decomposition: non-regular main diagonal block.");                              Esys_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getRILU:Break-down in RILU decomposition: non-regular main diagonal block.");
194                         }                         }
195                     }                     }
196                  }                  }
197                }                }
198             } /* end parallel region */             } /* end parallel region */
199    
200             if( Paso_noError()) {             if( Esys_noError()) {
201                /* if there are no nodes in the coarse level there is no more work to do */                /* if there are no nodes in the coarse level there is no more work to do */
202                out->n_C=n-out->n_F;                out->n_C=n-out->n_F;
203                if (out->n_C>0) {                if (out->n_C>0) {
204                     out->rows_in_C=MEMALLOC(out->n_C,index_t);                     out->rows_in_C=MEMALLOC(out->n_C,index_t);
205                     out->mask_C=MEMALLOC(n,index_t);                     out->mask_C=MEMALLOC(n,index_t);
206                     if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {                     if (! (Esys_checkPtr(out->mask_C) || Esys_checkPtr(out->rows_in_C) ) ) {
207                         /* creates an index for C from mask */                         /* creates an index for C from mask */
208                         #pragma omp parallel for private(i) schedule(static)                         #pragma omp parallel for private(i) schedule(static)
209                         for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];                         for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];
# Line 224  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 224  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
224                        } /* end parallel region */                        } /* end parallel region */
225                        /* get A_CF block: */                        /* get A_CF block: */
226                        out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F);                        out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F);
227                        if (Paso_noError()) {                        if (Esys_noError()) {
228                           /* get A_FC block: */                           /* get A_FC block: */
229                           out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);                           out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);
230                           /* get A_FF block: */                           /* get A_FF block: */
231                           if (Paso_noError()) {                           if (Esys_noError()) {
232                              schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C);                              schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C);
233                              time0=Paso_timer()-time0;                              time0=Esys_timer()-time0;
234                              if (Paso_noError()) {                              if (Esys_noError()) {
235                                  time1=Paso_timer();                                  time1=Esys_timer();
236                                  /* update A_CC block to get Schur complement and then apply RILU to it */                                  /* update A_CC block to get Schur complement and then apply RILU to it */
237                                  Paso_Solver_updateIncompleteSchurComplement(schur,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);                                  Paso_Solver_updateIncompleteSchurComplement(schur,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);
238                                  time1=Paso_timer()-time1;                                  time1=Esys_timer()-time1;
239                                  out->RILU_of_Schur=Paso_Solver_getRILU(schur,verbose);                                  out->RILU_of_Schur=Paso_Solver_getRILU(schur,verbose);
240                                  Paso_SparseMatrix_free(schur);                                  Paso_SparseMatrix_free(schur);
241                              }                              }
242                              /* allocate work arrays for RILU application */                              /* allocate work arrays for RILU application */
243                              if (Paso_noError()) {                              if (Esys_noError()) {
244                                out->x_F=MEMALLOC(n_block*out->n_F,double);                                out->x_F=MEMALLOC(n_block*out->n_F,double);
245                                out->b_F=MEMALLOC(n_block*out->n_F,double);                                out->b_F=MEMALLOC(n_block*out->n_F,double);
246                                out->x_C=MEMALLOC(n_block*out->n_C,double);                                out->x_C=MEMALLOC(n_block*out->n_C,double);
247                                out->b_C=MEMALLOC(n_block*out->n_C,double);                                out->b_C=MEMALLOC(n_block*out->n_C,double);
248                                if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {                                if (! (Esys_checkPtr(out->x_F) || Esys_checkPtr(out->b_F) || Esys_checkPtr(out->x_C) || Esys_checkPtr(out->b_C) ) ) {
249                                    #pragma omp parallel                                    #pragma omp parallel
250                                    {                                    {
251                                      #pragma omp for private(i,k) schedule(static)                                      #pragma omp for private(i,k) schedule(static)
# Line 275  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 275  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
275    }    }
276    TMPMEMFREE(mis_marker);    TMPMEMFREE(mis_marker);
277    TMPMEMFREE(counter);    TMPMEMFREE(counter);
278    if (Paso_noError()) {    if (Esys_noError()) {
279  /*  /*
280        if (verbose) {        if (verbose) {
281           printf("RILU: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F);           printf("RILU: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F);

Legend:
Removed from v.3232  
changed lines
  Added in v.3234

  ViewVC Help
Powered by ViewVC 1.1.26