← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to paso/src/AMLI.c

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
paso/src/AMLI.c
127
127
  
128
128
  /*Make sure we have block sizes 1*/
129
129
  if (A_p->col_block_size>1) {
130
 
     Paso_setError(TYPE_ERROR,"Paso_Solver_getAMLI: AMLI requires column block size 1.");
 
130
     Esys_setError(TYPE_ERROR,"Paso_Solver_getAMLI: AMLI requires column block size 1.");
131
131
     return NULL;
132
132
  }
133
133
  if (A_p->row_block_size>1) {
134
 
     Paso_setError(TYPE_ERROR,"Paso_Solver_getAMLI: AMLI requires row block size 1.");
 
134
     Esys_setError(TYPE_ERROR,"Paso_Solver_getAMLI: AMLI requires row block size 1.");
135
135
     return NULL;
136
136
  }
137
137
  out=MEMALLOC(1,Paso_Solver_AMLI);
138
138
  /* identify independend set of rows/columns */
139
139
  mis_marker=TMPMEMALLOC(n,index_t);
140
140
  counter=TMPMEMALLOC(n,index_t);
141
 
  if ( !( Paso_checkPtr(mis_marker) || Paso_checkPtr(counter) || Paso_checkPtr(out)) ) {
 
141
  if ( !( Esys_checkPtr(mis_marker) || Esys_checkPtr(counter) || Esys_checkPtr(out)) ) {
142
142
     out->AMLI_of_Schur=NULL;
143
143
     out->inv_A_FF=NULL;
144
144
     out->A_FF_pivot=NULL;
210
210
            }
211
211
        } else {
212
212
     
213
 
              if (Paso_noError()) {
 
213
              if (Esys_noError()) {
214
214
                 
215
215
                 /*#pragma omp parallel for private(i) schedule(static)
216
216
                 for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
228
228
                 out->rows_in_F=MEMALLOC(out->n_F,index_t);
229
229
                 out->inv_A_FF=MEMALLOC(n_block*n_block*out->n_F,double);
230
230
                 out->A_FF_pivot=NULL; /* later use for block size>3 */
231
 
                 if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->inv_A_FF) || Paso_checkPtr(out->rows_in_F) ) ) {
 
231
                 if (! (Esys_checkPtr(out->mask_F) || Esys_checkPtr(out->inv_A_FF) || Esys_checkPtr(out->rows_in_F) ) ) {
232
232
                    /* creates an index for F from mask */
233
233
                    #pragma omp parallel for private(i) schedule(static)
234
234
                    for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;
265
265
                    level=1;
266
266
               }
267
267
             
268
 
              if ( Paso_noError()) {
 
268
              if ( Esys_noError()) {
269
269
                    /* if there are no nodes in the coarse level there is no more work to do */
270
270
                    out->n_C=n-out->n_F;
271
271
                   
272
272
                   /*if (out->n_F>500) */
273
273
                    out->rows_in_C=MEMALLOC(out->n_C,index_t);
274
274
                    out->mask_C=MEMALLOC(n,index_t);
275
 
                    if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {
 
275
                    if (! (Esys_checkPtr(out->mask_C) || Esys_checkPtr(out->rows_in_C) ) ) {
276
276
                         /* creates an index for C from mask */
277
277
                         #pragma omp parallel for private(i) schedule(static)
278
278
                         for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];
290
290
                         }
291
291
                    }
292
292
              } 
293
 
              if ( Paso_noError()) {
 
293
              if ( Esys_noError()) {
294
294
                      /* get A_CF block: */
295
295
                      out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F);
296
296
                      /* get A_FC block: */
299
299
                      schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C);
300
300
      
301
301
              }
302
 
              if ( Paso_noError()) {
 
302
              if ( Esys_noError()) {
303
303
                     /*find the pattern of the schur complement with fill in*/
304
304
                    temp1=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern);
305
305
                    temp2=Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, temp1);
307
307
                    Paso_Pattern_free(temp1);
308
308
                    Paso_Pattern_free(temp2);
309
309
              }
310
 
              if ( Paso_noError()) {
 
310
              if ( Esys_noError()) {
311
311
                    /* copy values over*/ 
312
312
                    #pragma omp parallel for private(i,iPtr,j,iPtr_s,index,where_p) schedule(static)
313
313
                    for (i = 0; i < schur_withFillIn->numRows; ++i) {
329
329
                    out->AMLI_of_Schur=Paso_Solver_getAMLI(schur_withFillIn,level-1,options);
330
330
              }
331
331
              /* allocate work arrays for AMLI application */
332
 
              if (Paso_noError()) {
 
332
              if (Esys_noError()) {
333
333
                         out->x_F=MEMALLOC(n_block*out->n_F,double);
334
334
                         out->b_F=MEMALLOC(n_block*out->n_F,double);
335
335
                         out->x_C=MEMALLOC(n_block*out->n_C,double);
336
336
                         out->b_C=MEMALLOC(n_block*out->n_C,double);
337
337
      
338
 
                         if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {
 
338
                         if (! (Esys_checkPtr(out->x_F) || Esys_checkPtr(out->b_F) || Esys_checkPtr(out->x_C) || Esys_checkPtr(out->b_C) ) ) {
339
339
                             #pragma omp parallel for private(i) schedule(static)
340
340
                             for (i = 0; i < out->n_F; ++i) {
341
341
                                         out->x_F[i]=0.;
356
356
  TMPMEMFREE(mis_marker);
357
357
  TMPMEMFREE(counter);
358
358
 
359
 
  if (Paso_noError()) {
 
359
  if (Esys_noError()) {
360
360
      if (verbose && level>0 && !out->coarsest_level) {
361
361
         printf("AMLI: level: %d: %d unknowns eliminated. %d left.\n",level, out->n_F,out->n_C);
362
362
     }
410
410
     
411
411
     if (amli->coarsest_level) {
412
412
      
413
 
      time0=Paso_timer();
 
413
      time0=Esys_timer();
414
414
      /*If all unknown are eliminated then Jacobi is the best preconditioner*/
415
415
      if (amli->n_F==0 || amli->n_F==amli->n) {
416
416
         Paso_Preconditioner_LocalSmoother_solve(amli->A, amli->Smoother,x,b,1,FALSE);
428
428
         #endif
429
429
       #endif
430
430
       }
431
 
       time0=Paso_timer()-time0;
 
431
       time0=Esys_timer()-time0;
432
432
       if (verbose) fprintf(stderr,"timing: DIRECT SOLVER: %e\n",time0);
433
433
       
434
434
     } else {
435
435
        /* presmoothing */
436
 
         time0=Paso_timer();
 
436
         time0=Esys_timer();
437
437
         Paso_Preconditioner_LocalSmoother_solve(amli->A,amli->Smoother,x,b,pre_sweeps,FALSE);
438
 
         time0=Paso_timer()-time0;
 
438
         time0=Esys_timer()-time0;
439
439
         if (verbose) fprintf(stderr,"timing: Presmooting: %e\n",time0);
440
440
        /* end of presmoothing */
441
441
        
442
442
        
443
 
         time0=Paso_timer();
 
443
         time0=Esys_timer();
444
444
         #pragma omp parallel for private(i) schedule(static)
445
445
         for (i=0;i<amli->n;++i) r[i]=b[i];
446
446
         
461
461
        /* b_C=b_C-A_CF*x_F */
462
462
        Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amli->A_CF,amli->x_F,1.,amli->b_C);
463
463
        
464
 
        time0=Paso_timer()-time0;
 
464
        time0=Esys_timer()-time0;
465
465
        if (verbose) fprintf(stderr,"timing: Before next level: %e\n",time0);
466
466
        
467
467
        /* x_C=AMLI(b_C)     */
468
468
        Paso_Solver_solveAMLI(amli->AMLI_of_Schur,amli->x_C,amli->b_C);
469
469
        
470
 
        time0=Paso_timer();
 
470
        time0=Esys_timer();
471
471
        
472
472
        /* b_F=-A_FC*x_C */ 
473
473
        Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amli->A_FC,amli->x_C,0.,amli->b_F);
490
490
            }
491
491
        }
492
492
        
493
 
        time0=Paso_timer()-time0;
 
493
        time0=Esys_timer()-time0;
494
494
        if (verbose) fprintf(stderr,"timing: After next level: %e\n",time0);
495
495
 
496
496
     /*postsmoothing*/
497
 
     time0=Paso_timer();
 
497
     time0=Esys_timer();
498
498
     Paso_Preconditioner_LocalSmoother_solve(amli->A,amli->Smoother,x,b,post_sweeps,TRUE);     
499
 
     time0=Paso_timer()-time0;
 
499
     time0=Esys_timer()-time0;
500
500
     if (verbose) fprintf(stderr,"timing: Postsmoothing: %e\n",time0);
501
501
 
502
502
     /*end of postsmoothing*/