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

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

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

revision 1638 by phornby, Fri Jul 11 13:12:46 2008 UTC revision 1639 by gross, Mon Jul 14 08:55:25 2008 UTC
# Line 57  void Paso_Solver(Paso_SystemMatrix* A,do Line 57  void Paso_Solver(Paso_SystemMatrix* A,do
57     dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);     dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
58     double blocktimer_start = blocktimer_time();     double blocktimer_start = blocktimer_time();
59    
60     tolerance=MAX(options->tolerance,EPSILON);       tolerance=MAX(options->tolerance,EPSILON);
61     Paso_resetError();       Paso_resetError();
62     method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);       method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
63     /* check matrix type */       /* check matrix type */
64     if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {       if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
65        Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");         Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
66     }       }
67     if (A->col_block_size != A->row_block_size) {       if (A->col_block_size != A->row_block_size) {
68        Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal.");          Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal.");
69     }       }
70     if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {       if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {
71        Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");          Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
72        return;          return;
73     }       }
74     Performance_startMonitor(pp,PERFORMANCE_ALL);       /* this for testing only */
75     if (Paso_noError()) {       if (method==PASO_NONLINEAR_GMRES) {
76        /* get normalization */          Paso_Function* F=NULL;
77        scaling=Paso_SystemMatrix_borrowNormalization(A);          F=Paso_Function_LinearSystem_alloc(A,b,options);
78        if (Paso_noError()) {          errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);
79           /* get the norm of the right hand side */          if (errorCode==NO_ERROR) {
80           norm2_of_b=0.;             Paso_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");
81           norm_max_of_b=0.;          }
82  #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)          Paso_Function_LinearSystem_free(F);
83           {       }
84              norm2_of_b_local=0.;       /* ========================= */
85              norm_max_of_b_local=0.;       Performance_startMonitor(pp,PERFORMANCE_ALL);
86  #pragma omp for private(i) schedule(static)       if (Paso_noError()) {
87              for (i = 0; i < numEqua ; ++i) {          /* get normalization */
88                 norm2_of_b_local += b[i] * b[i];          scaling=Paso_SystemMatrix_borrowNormalization(A);
89                 norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);          if (Paso_noError()) {
90              }             /* get the norm of the right hand side */
91  #pragma omp critical             norm2_of_b=0.;
92              {             norm_max_of_b=0.;
93                 norm2_of_b += norm2_of_b_local;             #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
94                 norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);             {
95              }                  norm2_of_b_local=0.;
96           }                  norm_max_of_b_local=0.;
97           /* TODO: use one call */                  #pragma omp for private(i) schedule(static)
98  #ifdef PASO_MPI                  for (i = 0; i < numEqua ; ++i) {
99           {                        norm2_of_b_local += b[i] * b[i];
100              loc_norm = norm2_of_b;                        norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
101              MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);                  }
102              loc_norm = norm_max_of_b;                  #pragma omp critical
103              MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);                  {
104           }                     norm2_of_b += norm2_of_b_local;
105  #endif                     norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
106           norm2_of_b=sqrt(norm2_of_b);                  }
107               }
108               /* TODO: use one call */
109               #ifdef PASO_MPI
110               {
111                   loc_norm = norm2_of_b;
112                   MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
113                   loc_norm = norm_max_of_b;
114                   MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
115               }
116               #endif
117               norm2_of_b=sqrt(norm2_of_b);
118            
119           /* if norm2_of_b==0 we are ready: x=0 */           /* if norm2_of_b==0 we are ready: x=0 */
120           if (norm2_of_b <=0.) {           if (norm2_of_b <=0.) {
# Line 239  void Paso_Solver(Paso_SystemMatrix* A,do Line 250  void Paso_Solver(Paso_SystemMatrix* A,do
250                MEMFREE(r);                MEMFREE(r);
251                time_iter=Paso_timer()-time_iter;                time_iter=Paso_timer()-time_iter;
252                if (options->verbose)  {                if (options->verbose)  {
253                   printf("timing: solver: %.4e sec\n",time_iter);                   printf("\ntiming: solver: %.4e sec\n",time_iter);
254                   if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);                   if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
255                }                }
256             }             }

Legend:
Removed from v.1638  
changed lines
  Added in v.1639

  ViewVC Help
Powered by ViewVC 1.1.26