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

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

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

revision 1872 by jfenwick, Mon Oct 13 00:18:55 2008 UTC revision 1902 by artak, Wed Oct 22 03:54:14 2008 UTC
# Line 408  void Paso_Solver_solveGS(Paso_Solver_GS Line 408  void Paso_Solver_solveGS(Paso_Solver_GS
408       return;       return;
409  }  }
410    
411    /**************************************************************/
412    
413    /* apply GS precondition b-> x                              
414    
415         in fact it solves LD^{-1}Ux=b in the form x= U^{-1} D L^{-1}b
416    
417     should be called within a parallel region                                              
418     barrier synconization should be performed to make sure that the input vector available
419    
420    */
421    
422    void Paso_Solver_solveGS1(Paso_Solver_GS * gs, double * x, double * b) {
423         register dim_t i,k;
424         register index_t color,iptr_ik,iptr_main;
425         register double A11,A12,A21,A22,A13,A23,A33,A32,A31,S1,S2,S3,R1,R2,R3,D,S21,S22,S12,S11,S13,S23,S33,S32,S31;
426         dim_t n_block=gs->n_block;
427         dim_t n=gs->n;
428         index_t* pivot=NULL;
429        
430         /* copy x into b*/
431         #pragma omp parallel for private(i) schedule(static)
432         /*for (i=0;i<n*n_block;++i) x[i]=b[i];*/
433        
434         /* forward substitution */
435         for (color=0;color<gs->num_colors;++color) {
436               if (n_block==1) {
437                  #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1,iptr_main)
438                  for (i = 0; i < n; ++i) {
439                       if (gs->colorOf[i]==color) {
440                         /* x_i=x_i-a_ik*x_k */                    
441                         S1=b[i];
442                         for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {
443                              k=gs->pattern->index[iptr_ik];                          
444                              if (gs->colorOf[k]<color) {
445                                 R1=x[k];                              
446                                 S1-=gs->factors->val[iptr_ik]*R1;
447                               }
448                         }
449                         iptr_main=gs->main_iptr[i];
450                         x[i]=(1/gs->factors->val[iptr_main])*S1;
451                       }
452                  }
453               } else if (n_block==2) {
454                  #pragma omp parallel for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)
455                  for (i = 0; i < n; ++i) {
456                       if (gs->colorOf[i]==color) {
457                         /* x_i=x_i-a_ik*x_k */
458                         S1=b[2*i];
459                         S2=b[2*i+1];
460                         for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {
461                              k=gs->pattern->index[iptr_ik];                          
462                              if (gs->colorOf[k]<color) {
463                                 R1=x[2*k];
464                                 R2=x[2*k+1];
465                                 S1-=gs->factors->val[4*iptr_ik  ]*R1+gs->factors->val[4*iptr_ik+2]*R2;
466                                 S2-=gs->factors->val[4*iptr_ik+1]*R1+gs->factors->val[4*iptr_ik+3]*R2;
467                              }
468                         }
469                         iptr_main=gs->main_iptr[i];
470                         A11=gs->factors->val[iptr_main*4];
471                         A21=gs->factors->val[iptr_main*4+1];
472                         A12=gs->factors->val[iptr_main*4+2];
473                         A22=gs->factors->val[iptr_main*4+3];
474                         D = A11*A22-A12*A21;
475                         if (ABS(D)>0.) {
476                              D=1./D;
477                              S11= A22*D;
478                              S21=-A21*D;
479                              S12=-A12*D;
480                              S22= A11*D;
481                              x[2*i  ]=S11*S1+S12*S2;
482                              x[2*i+1]=S21*S1+S22*S2;
483                         } else {
484                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");
485                           }
486                       }
487    
488                  }
489               } else if (n_block==3) {
490                  #pragma omp parallel for schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3)
491                  for (i = 0; i < n; ++i) {
492                       if (gs->colorOf[i]==color) {
493                         /* x_i=x_i-a_ik*x_k */
494                         S1=b[3*i];
495                         S2=b[3*i+1];
496                         S3=b[3*i+2];
497                         for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {
498                              k=gs->pattern->index[iptr_ik];                          
499                              if (gs->colorOf[k]<color) {
500                                 R1=x[3*k];
501                                 R2=x[3*k+1];
502                                 R3=x[3*k+2];
503                                 S1-=gs->factors->val[9*iptr_ik  ]*R1+gs->factors->val[9*iptr_ik+3]*R2+gs->factors->val[9*iptr_ik+6]*R3;
504                                 S2-=gs->factors->val[9*iptr_ik+1]*R1+gs->factors->val[9*iptr_ik+4]*R2+gs->factors->val[9*iptr_ik+7]*R3;
505                                 S3-=gs->factors->val[9*iptr_ik+2]*R1+gs->factors->val[9*iptr_ik+5]*R2+gs->factors->val[9*iptr_ik+8]*R3;
506                              }
507                         }
508                         iptr_main=gs->main_iptr[i];
509                         A11=gs->factors->val[iptr_main*9  ];
510                         A21=gs->factors->val[iptr_main*9+1];
511                         A31=gs->factors->val[iptr_main*9+2];
512                         A12=gs->factors->val[iptr_main*9+3];
513                         A22=gs->factors->val[iptr_main*9+4];
514                         A32=gs->factors->val[iptr_main*9+5];
515                         A13=gs->factors->val[iptr_main*9+6];
516                         A23=gs->factors->val[iptr_main*9+7];
517                         A33=gs->factors->val[iptr_main*9+8];
518                         D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
519                         if (ABS(D)>0.) {
520                              D=1./D;
521                              S11=(A22*A33-A23*A32)*D;
522                              S21=(A31*A23-A21*A33)*D;
523                              S31=(A21*A32-A31*A22)*D;
524                              S12=(A13*A32-A12*A33)*D;
525                              S22=(A11*A33-A31*A13)*D;
526                              S32=(A12*A31-A11*A32)*D;
527                              S13=(A12*A23-A13*A22)*D;
528                              S23=(A13*A21-A11*A23)*D;
529                              S33=(A11*A22-A12*A21)*D;
530                                 /* a_ik=a_ii^{-1}*a_ik */
531                              x[3*i  ]=S11*S1+S12*S2+S13*S3;
532                              x[3*i+1]=S21*S1+S22*S2+S23*S3;
533                              x[3*i+2]=S31*S1+S32*S2+S33*S3;
534                           } else {
535                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");
536                           }
537                    }
538                  }
539               }
540         }
541        
542         /* Multipling with diag(A) */
543         Paso_Solver_applyBlockDiagonalMatrix(gs->n_block,gs->n,gs->diag,pivot,x,x);
544    
545         /* backward substitution */
546         for (color=(gs->num_colors)-1;color>-1;--color) {
547               if (n_block==1) {
548                  #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1)
549                  for (i = 0; i < n; ++i) {
550                       if (gs->colorOf[i]==color) {
551                         /* x_i=x_i-a_ik*x_k */
552                         S1=x[i];
553                         for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {
554                              k=gs->pattern->index[iptr_ik];                          
555                              if (gs->colorOf[k]>color) {
556                                 R1=x[k];
557                                 S1-=gs->factors->val[iptr_ik]*R1;
558                              }
559                         }
560                         /*x[i]=S1;*/
561                         iptr_main=gs->main_iptr[i];
562                         x[i]=(1/gs->factors->val[iptr_main])*S1;
563                       }
564                  }
565               } else if (n_block==2) {
566                  #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)
567                  for (i = 0; i < n; ++i) {
568                       if (gs->colorOf[i]==color) {
569                         /* x_i=x_i-a_ik*x_k */
570                         S1=x[2*i];
571                         S2=x[2*i+1];
572                         for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {
573                              k=gs->pattern->index[iptr_ik];                          
574                              if (gs->colorOf[k]>color) {
575                                 R1=x[2*k];
576                                 R2=x[2*k+1];
577                                 S1-=gs->factors->val[4*iptr_ik  ]*R1+gs->factors->val[4*iptr_ik+2]*R2;
578                                 S2-=gs->factors->val[4*iptr_ik+1]*R1+gs->factors->val[4*iptr_ik+3]*R2;
579                              }
580                         }
581                         /*x[2*i]=S1;
582                         x[2*i+1]=S2;*/
583                         iptr_main=gs->main_iptr[i];
584                         A11=gs->factors->val[iptr_main*4];
585                         A21=gs->factors->val[iptr_main*4+1];
586                         A12=gs->factors->val[iptr_main*4+2];
587                         A22=gs->factors->val[iptr_main*4+3];
588                         D = A11*A22-A12*A21;
589                         if (ABS(D)>0.) {
590                              D=1./D;
591                              S11= A22*D;
592                              S21=-A21*D;
593                              S12=-A12*D;
594                              S22= A11*D;
595                              x[2*i  ]=S11*S1+S12*S2;
596                              x[2*i+1]=S21*S1+S22*S2;
597                         } else {
598                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");
599                           }
600    
601                        }
602                  }
603               } else if (n_block==3) {
604                  #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)
605                  for (i = 0; i < n; ++i) {
606                       if (gs->colorOf[i]==color) {
607                         /* x_i=x_i-a_ik*x_k */
608                         S1=x[3*i  ];
609                         S2=x[3*i+1];
610                         S3=x[3*i+2];
611                         for (iptr_ik=gs->pattern->ptr[i];iptr_ik<gs->pattern->ptr[i+1]; ++iptr_ik) {
612                              k=gs->pattern->index[iptr_ik];                          
613                              if (gs->colorOf[k]>color) {
614                                 R1=x[3*k];
615                                 R2=x[3*k+1];
616                                 R3=x[3*k+2];
617                                 S1-=gs->factors->val[9*iptr_ik  ]*R1+gs->factors->val[9*iptr_ik+3]*R2+gs->factors->val[9*iptr_ik+6]*R3;
618                                 S2-=gs->factors->val[9*iptr_ik+1]*R1+gs->factors->val[9*iptr_ik+4]*R2+gs->factors->val[9*iptr_ik+7]*R3;
619                                 S3-=gs->factors->val[9*iptr_ik+2]*R1+gs->factors->val[9*iptr_ik+5]*R2+gs->factors->val[9*iptr_ik+8]*R3;
620                              }
621                         }
622    /*                     x[3*i]=S1;
623                         x[3*i+1]=S2;
624                         x[3*i+2]=S3;
625    */                   iptr_main=gs->main_iptr[i];
626                         A11=gs->factors->val[iptr_main*9  ];
627                         A21=gs->factors->val[iptr_main*9+1];
628                         A31=gs->factors->val[iptr_main*9+2];
629                         A12=gs->factors->val[iptr_main*9+3];
630                         A22=gs->factors->val[iptr_main*9+4];
631                         A32=gs->factors->val[iptr_main*9+5];
632                         A13=gs->factors->val[iptr_main*9+6];
633                         A23=gs->factors->val[iptr_main*9+7];
634                         A33=gs->factors->val[iptr_main*9+8];
635                         D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
636                         if (ABS(D)>0.) {
637                              D=1./D;
638                              S11=(A22*A33-A23*A32)*D;
639                              S21=(A31*A23-A21*A33)*D;
640                              S31=(A21*A32-A31*A22)*D;
641                              S12=(A13*A32-A12*A33)*D;
642                              S22=(A11*A33-A31*A13)*D;
643                              S32=(A12*A31-A11*A32)*D;
644                              S13=(A12*A23-A13*A22)*D;
645                              S23=(A13*A21-A11*A23)*D;
646                              S33=(A11*A22-A12*A21)*D;
647                              x[3*i  ]=S11*S1+S12*S2+S13*S3;
648                              x[3*i+1]=S21*S1+S22*S2+S23*S3;
649                              x[3*i+2]=S31*S1+S32*S2+S33*S3;
650                           } else {
651                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");
652                           }
653                       }
654                  }
655             }
656         }
657        
658         return;
659    }
660    

Legend:
Removed from v.1872  
changed lines
  Added in v.1902

  ViewVC Help
Powered by ViewVC 1.1.26