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

Revision 1787 - (hide annotations)
Mon Sep 15 01:36:34 2008 UTC (10 years, 7 months ago) by artak
File MIME type: text/plain
File size: 12961 byte(s)
```MINRES solver is added to escript. Additional 16 tests are added to run_simplesolve for MINRES and TFQMR solvers
```
 1 ksteube 1312 2 jgs 150 /* \$Id\$ */ 3 4 ksteube 1312 /******************************************************* 5 * 6 * Copyright 2003-2007 by ACceSS MNRF 7 * Copyright 2007 by University of Queensland 8 * 9 * http://esscc.uq.edu.au 10 * Primary Business: Queensland, Australia 11 * Licensed under the Open Software License version 3.0 12 * http://www.opensource.org/licenses/osl-3.0.php 13 * 14 *******************************************************/ 15 dhawcroft 631 16 jgs 150 /**************************************************************/ 17 18 /* Paso: SystemMatrix: controls iterative linear system solvers */ 19 20 /**************************************************************/ 21 22 /* Copyrights by ACcESS Australia 2003/04 */ 23 /* Author: gross@access.edu.au */ 24 25 /**************************************************************/ 26 27 gross 700 #include "Paso.h" 28 #include "SystemMatrix.h" 29 jgs 150 #include "Solver.h" 30 ksteube 1312 #include "escript/blocktimer.h" 31 jgs 150 32 /***********************************************************************************/ 33 34 /* free space */ 35 36 void Paso_Solver_free(Paso_SystemMatrix* A) { 37 phornby 1628 Paso_Preconditioner_free(A->solver); 38 A->solver=NULL; 39 jgs 150 } 40 /* call the iterative solver: */ 41 42 phornby 1628 void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b, 43 Paso_Options* options,Paso_Performance* pp) { 44 45 double norm2_of_b,tol,tolerance,time_iter; 46 double *r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b; 47 double norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local; 48 double norm_max_of_residual_local,norm_max_of_residual; 49 double last_norm_max_of_residual,*scaling; 50 #ifdef PASO_MPI 51 double loc_norm; 52 #endif 53 dim_t i,totIter,cntIter,method; 54 bool_t finalizeIteration; 55 err_t errorCode; 56 dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A); 57 dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A); 58 ksteube 1707 double blocktimer_precond, blocktimer_start = blocktimer_time(); 59 gross 415 60 gross 1639 tolerance=MAX(options->tolerance,EPSILON); 61 Paso_resetError(); 62 method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric); 63 /* check matrix type */ 64 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."); 66 } 67 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."); 69 } 70 if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) { 71 artak 1787 Paso_setError(TYPE_ERROR,"Iterative solver requires a square matrix."); 72 gross 1639 return; 73 } 74 /* this for testing only */ 75 if (method==PASO_NONLINEAR_GMRES) { 76 Paso_Function* F=NULL; 77 F=Paso_Function_LinearSystem_alloc(A,b,options); 78 errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp); 79 if (errorCode==NO_ERROR) { 80 Paso_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured."); 81 } 82 Paso_Function_LinearSystem_free(F); 83 } 84 /* ========================= */ 85 Performance_startMonitor(pp,PERFORMANCE_ALL); 86 if (Paso_noError()) { 87 /* get normalization */ 88 scaling=Paso_SystemMatrix_borrowNormalization(A); 89 if (Paso_noError()) { 90 /* get the norm of the right hand side */ 91 norm2_of_b=0.; 92 norm_max_of_b=0.; 93 #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local) 94 { 95 norm2_of_b_local=0.; 96 norm_max_of_b_local=0.; 97 #pragma omp for private(i) schedule(static) 98 for (i = 0; i < numEqua ; ++i) { 99 norm2_of_b_local += b[i] * b[i]; 100 norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local); 101 } 102 #pragma omp critical 103 { 104 norm2_of_b += norm2_of_b_local; 105 norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b); 106 } 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 jgs 150 119 phornby 1628 /* if norm2_of_b==0 we are ready: x=0 */ 120 if (norm2_of_b <=0.) { 121 #pragma omp parallel for private(i) schedule(static) 122 for (i = 0; i < numSol; i++) x[i]=0.; 123 if (options->verbose) printf("right hand side is identical zero.\n"); 124 } else { 125 if (options->verbose) { 126 printf("l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b); 127 printf("l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance); 128 switch (method) { 129 case PASO_BICGSTAB: 130 printf("Iterative method is BiCGStab.\n"); 131 break; 132 case PASO_PCG: 133 printf("Iterative method is PCG.\n"); 134 break; 135 artak 1703 case PASO_TFQMR: 136 printf("Iterative method is TFQMR.\n"); 137 break; 138 artak 1787 case PASO_MINRES: 139 printf("Iterative method is MINRES.\n"); 140 break; 141 phornby 1628 case PASO_PRES20: 142 printf("Iterative method is PRES20.\n"); 143 break; 144 case PASO_GMRES: 145 if (options->restart>0) { 146 printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart); 147 } else { 148 printf("Iterative method is GMRES(%d)\n",options->truncation); 149 gross 584 } 150 phornby 1628 break; 151 } 152 } 153 /* construct the preconditioner */ 154 gross 584 155 ksteube 1707 blocktimer_precond = blocktimer_time(); 156 phornby 1628 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT); 157 Paso_Solver_setPreconditioner(A,options); 158 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT); 159 ksteube 1707 blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond); 160 phornby 1628 if (! Paso_noError()) return; 161 ksteube 1312 time_iter=Paso_timer(); 162 gross 584 /* get an initial guess by evaluating the preconditioner */ 163 gross 1556 Paso_Solver_solvePreconditioner(A,x,b); 164 gross 584 /* start the iteration process :*/ 165 ksteube 1312 r=TMPMEMALLOC(numEqua,double); 166 Paso_checkPtr(r); 167 if (Paso_noError()) { 168 totIter = 0; 169 finalizeIteration = FALSE; 170 last_norm2_of_residual=norm2_of_b; 171 last_norm_max_of_residual=norm_max_of_b; 172 artak 1787 /* Loop */ 173 ksteube 1312 while (! finalizeIteration) { 174 cntIter = options->iter_max - totIter; 175 finalizeIteration = TRUE; 176 /* Set initial residual. */ 177 norm2_of_residual = 0; 178 norm_max_of_residual = 0; 179 gross 1556 #pragma omp parallel for private(i) schedule(static) 180 for (i = 0; i < numEqua; i++) r[i]=b[i]; 181 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r); 182 183 ksteube 1312 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local) 184 { 185 norm2_of_residual_local = 0; 186 norm_max_of_residual_local = 0; 187 #pragma omp for private(i) schedule(static) 188 for (i = 0; i < numEqua; i++) { 189 norm2_of_residual_local+= r[i] * r[i]; 190 norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local); 191 } 192 #pragma omp critical 193 { 194 norm2_of_residual += norm2_of_residual_local; 195 norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual); 196 } 197 gross 584 } 198 ksteube 1312 /* TODO: use one call */ 199 #ifdef PASO_MPI 200 loc_norm = norm2_of_residual; 201 MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm); 202 loc_norm = norm_max_of_residual; 203 MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm); 204 #endif 205 norm2_of_residual =sqrt(norm2_of_residual); 206 207 phornby 1628 if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual); 208 if (totIter>0 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) { 209 if (options->verbose) printf(" divergence!\n"); 210 Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up."); 211 } else { 212 /* */ 213 if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) { 214 tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b); 215 if (options->verbose) printf(" (new tolerance = %e).\n",tol); 216 last_norm2_of_residual=norm2_of_residual; 217 last_norm_max_of_residual=norm_max_of_residual; 218 /* call the solver */ 219 switch (method) { 220 case PASO_BICGSTAB: 221 errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp); 222 break; 223 case PASO_PCG: 224 errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp); 225 break; 226 artak 1703 case PASO_TFQMR: 227 artak 1787 errorCode = Paso_Solver_TFQMR(A, r, x, &cntIter, &tol, pp); 228 artak 1703 break; 229 artak 1787 case PASO_MINRES: 230 errorCode = Paso_Solver_MINRES(A, r, x, &cntIter, &tolerance, &tol, pp); 231 break; 232 phornby 1628 case PASO_PRES20: 233 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp); 234 break; 235 case PASO_GMRES: 236 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp); 237 break; 238 } 239 totIter += cntIter; 240 /* error handling */ 241 if (errorCode==NO_ERROR) { 242 finalizeIteration = FALSE; 243 } else if (errorCode==SOLVER_MAXITER_REACHED) { 244 Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion."); 245 if (options->verbose) printf("Maximum number of iterations reached.!\n"); 246 } else if (errorCode == SOLVER_INPUT_ERROR ) { 247 Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver."); 248 if (options->verbose) printf("Internal error!\n"); 249 } else if ( errorCode == SOLVER_BREAKDOWN ) { 250 if (cntIter <= 1) { 251 Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver."); 252 if (options->verbose) printf("Uncurable break down!\n"); 253 } else { 254 if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol); 255 finalizeIteration = FALSE; 256 ksteube 1312 } 257 } else { 258 if (options->verbose) printf(". convergence! \n"); 259 } 260 jgs 150 } 261 ksteube 1312 } /* while */ 262 } 263 gross 584 MEMFREE(r); 264 time_iter=Paso_timer()-time_iter; 265 if (options->verbose) { 266 gross 1639 printf("\ntiming: solver: %.4e sec\n",time_iter); 267 gross 584 if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter); 268 jgs 150 } 269 } 270 } 271 gross 584 } 272 phornby 1628 } 273 Performance_stopMonitor(pp,PERFORMANCE_ALL); 274 blocktimer_increment("Paso_Solver()", blocktimer_start); 275 jgs 150 }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision