# Contents of /trunk/paso/src/Solver.c

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

## Properties

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