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

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

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

revision 1822 by artak, Tue Sep 30 05:58:06 2008 UTC revision 1823 by artak, Wed Oct 1 05:56:05 2008 UTC
# Line 67  Paso_Solver_GS* Paso_Solver_getGS(Paso_S Line 67  Paso_Solver_GS* Paso_Solver_getGS(Paso_S
67    out->factors=Paso_SparseMatrix_getReference(A);    out->factors=Paso_SparseMatrix_getReference(A);
68    out->n_block=n_block;    out->n_block=n_block;
69    out->n=n;    out->n=n;
70      
71    if ( !(Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->factors)) ) {    if ( !(Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->factors)) ) {
72      time0=Paso_timer();      time0=Paso_timer();
73      Paso_Pattern_color(A->pattern,&out->num_colors,out->colorOf);      Paso_Pattern_color(A->pattern,&out->num_colors,out->colorOf);
# Line 173  Paso_Solver_GS* Paso_Solver_getGS(Paso_S Line 173  Paso_Solver_GS* Paso_Solver_getGS(Paso_S
173  void Paso_Solver_solveGS(Paso_Solver_GS * gs, double * x, double * b) {  void Paso_Solver_solveGS(Paso_Solver_GS * gs, double * x, double * b) {
174       register dim_t i,k;       register dim_t i,k;
175       register index_t color,iptr_ik,iptr_main;       register index_t color,iptr_ik,iptr_main;
176       register double S1,S2,S3,R1,R2,R3;       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;
177       dim_t n_block=gs->n_block;       dim_t n_block=gs->n_block;
178       dim_t n=gs->n;       dim_t n=gs->n;
179       index_t* pivot=NULL;       index_t* pivot=NULL;
# Line 217  void Paso_Solver_solveGS(Paso_Solver_GS Line 217  void Paso_Solver_solveGS(Paso_Solver_GS
217                            }                            }
218                       }                       }
219                       iptr_main=gs->main_iptr[i];                       iptr_main=gs->main_iptr[i];
220                       x[2*i  ]=(1/gs->factors->val[4*iptr_main  ])*S1+(1/gs->factors->val[4*iptr_main+2])*S2;                       A11=gs->factors->val[iptr_main*4];
221                       x[2*i+1]=(1/gs->factors->val[4*iptr_main+1])*S1+(1/gs->factors->val[4*iptr_main+3])*S2;                       A21=gs->factors->val[iptr_main*4+1];
222                         A12=gs->factors->val[iptr_main*4+2];
223                         A22=gs->factors->val[iptr_main*4+3];
224                         D = A11*A22-A12*A21;
225                         if (ABS(D)>0.) {
226                              D=1./D;
227                              S11= A22*D;
228                              S21=-A21*D;
229                              S12=-A12*D;
230                              S22= A11*D;
231                              x[2*i  ]=S11*S1+S12*S2;
232                              x[2*i+1]=S21*S1+S22*S2;
233                         } else {
234                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");
235                           }
236                     }                     }
237    
238                }                }
# Line 242  void Paso_Solver_solveGS(Paso_Solver_GS Line 256  void Paso_Solver_solveGS(Paso_Solver_GS
256                            }                            }
257                       }                       }
258                       iptr_main=gs->main_iptr[i];                       iptr_main=gs->main_iptr[i];
259                       x[3*i  ]=(1/gs->factors->val[9*iptr_main  ])*S1+(1/gs->factors->val[9*iptr_main+3])*S2+(1/gs->factors->val[9*iptr_main+6])*S3;                       A11=gs->factors->val[iptr_main*9  ];
260                       x[3*i+1]=(1/gs->factors->val[9*iptr_main+1])*S1+(1/gs->factors->val[9*iptr_main+4])*S2+(1/gs->factors->val[9*iptr_main+7])*S3;                       A21=gs->factors->val[iptr_main*9+1];
261                       x[3*i+2]=(1/gs->factors->val[9*iptr_main+2])*S1+(1/gs->factors->val[9*iptr_main+5])*S2+(1/gs->factors->val[9*iptr_main+8])*S3;                       A31=gs->factors->val[iptr_main*9+2];
262                   }                       A12=gs->factors->val[iptr_main*9+3];
263                         A22=gs->factors->val[iptr_main*9+4];
264                         A32=gs->factors->val[iptr_main*9+5];
265                         A13=gs->factors->val[iptr_main*9+6];
266                         A23=gs->factors->val[iptr_main*9+7];
267                         A33=gs->factors->val[iptr_main*9+8];
268                         D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
269                         if (ABS(D)>0.) {
270                              D=1./D;
271                              S11=(A22*A33-A23*A32)*D;
272                              S21=(A31*A23-A21*A33)*D;
273                              S31=(A21*A32-A31*A22)*D;
274                              S12=(A13*A32-A12*A33)*D;
275                              S22=(A11*A33-A31*A13)*D;
276                              S32=(A12*A31-A11*A32)*D;
277                              S13=(A12*A23-A13*A22)*D;
278                              S23=(A13*A21-A11*A23)*D;
279                              S33=(A11*A22-A12*A21)*D;
280                                 /* a_ik=a_ii^{-1}*a_ik */
281                              x[3*i  ]=S11*S1+S12*S2+S13*S3;
282                              x[3*i+1]=S21*S1+S22*S2+S23*S3;
283                              x[3*i+2]=S31*S1+S32*S2+S33*S3;
284                           } else {
285                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");
286                           }
287                    }
288                }                }
289             }             }
290       }       }
# Line 292  void Paso_Solver_solveGS(Paso_Solver_GS Line 331  void Paso_Solver_solveGS(Paso_Solver_GS
331                       /*x[2*i]=S1;                       /*x[2*i]=S1;
332                       x[2*i+1]=S2;*/                       x[2*i+1]=S2;*/
333                       iptr_main=gs->main_iptr[i];                       iptr_main=gs->main_iptr[i];
334                       x[2*i  ]=(1/gs->factors->val[4*iptr_main  ])*S1+(1/gs->factors->val[4*iptr_main+2])*S2;                       A11=gs->factors->val[iptr_main*4];
335                       x[2*i+1]=(1/gs->factors->val[4*iptr_main+1])*S1+(1/gs->factors->val[4*iptr_main+3])*S2;                       A21=gs->factors->val[iptr_main*4+1];
336                     }                       A12=gs->factors->val[iptr_main*4+2];
337                         A22=gs->factors->val[iptr_main*4+3];
338                         D = A11*A22-A12*A21;
339                         if (ABS(D)>0.) {
340                              D=1./D;
341                              S11= A22*D;
342                              S21=-A21*D;
343                              S12=-A12*D;
344                              S22= A11*D;
345                              x[2*i  ]=S11*S1+S12*S2;
346                              x[2*i+1]=S21*S1+S22*S2;
347                         } else {
348                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");
349                           }
350    
351                        }
352                }                }
353             } else if (n_block==3) {             } else if (n_block==3) {
354                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)
# Line 318  void Paso_Solver_solveGS(Paso_Solver_GS Line 372  void Paso_Solver_solveGS(Paso_Solver_GS
372  /*                     x[3*i]=S1;  /*                     x[3*i]=S1;
373                       x[3*i+1]=S2;                       x[3*i+1]=S2;
374                       x[3*i+2]=S3;                       x[3*i+2]=S3;
375  */                     iptr_main=gs->main_iptr[i];  */                   iptr_main=gs->main_iptr[i];
376                       x[3*i  ]=(1/gs->factors->val[9*iptr_main  ])*S1+(1/gs->factors->val[9*iptr_main+3])*S2+(1/gs->factors->val[9*iptr_main+6])*S3;                       A11=gs->factors->val[iptr_main*9  ];
377                       x[3*i+1]=(1/gs->factors->val[9*iptr_main+1])*S1+(1/gs->factors->val[9*iptr_main+4])*S2+(1/gs->factors->val[9*iptr_main+7])*S3;                       A21=gs->factors->val[iptr_main*9+1];
378                       x[3*i+2]=(1/gs->factors->val[9*iptr_main+2])*S1+(1/gs->factors->val[9*iptr_main+5])*S2+(1/gs->factors->val[9*iptr_main+8])*S3;                       A31=gs->factors->val[iptr_main*9+2];
379                                             A12=gs->factors->val[iptr_main*9+3];
380                         A22=gs->factors->val[iptr_main*9+4];
381                         A32=gs->factors->val[iptr_main*9+5];
382                         A13=gs->factors->val[iptr_main*9+6];
383                         A23=gs->factors->val[iptr_main*9+7];
384                         A33=gs->factors->val[iptr_main*9+8];
385                         D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
386                         if (ABS(D)>0.) {
387                              D=1./D;
388                              S11=(A22*A33-A23*A32)*D;
389                              S21=(A31*A23-A21*A33)*D;
390                              S31=(A21*A32-A31*A22)*D;
391                              S12=(A13*A32-A12*A33)*D;
392                              S22=(A11*A33-A31*A13)*D;
393                              S32=(A12*A31-A11*A32)*D;
394                              S13=(A12*A23-A13*A22)*D;
395                              S23=(A13*A21-A11*A23)*D;
396                              S33=(A11*A22-A12*A21)*D;
397                              x[3*i  ]=S11*S1+S12*S2+S13*S3;
398                              x[3*i+1]=S21*S1+S22*S2+S23*S3;
399                              x[3*i+2]=S31*S1+S32*S2+S33*S3;
400                           } else {
401                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getGS: non-regular main diagonal block.");
402                           }
403                     }                     }
404                }                }
405           }           }
406       }       }
407        
408         if (gs->sweeps>1) {
409         /* Compute the residual b=b-Ax*/
410         Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), gs->factors, x, DBLE(2), b);
411         /* Go round again*/
412         gs->sweeps=gs->sweeps-1;
413         Paso_Solver_solveGS(gs,x,b);
414         }
415        
416       return;       return;
417  }  }
418    
Legend:
Removed from v.1822  
changed lines
  Added in v.1823
  ViewVC Help
Powered by ViewVC 1.1.26