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

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

## Properties

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