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

Revision 1384 - (show annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 3 months ago) by phornby
Original Path: temp_trunk_copy/paso/src/Solver.c
File MIME type: text/plain
File size: 12938 byte(s)
Make a temp copy of the trunk before checking in the windows changes

 1 2 /* \$Id\$ */ 3 4 /******************************************************* 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 16 /**************************************************************/ 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 #include "Paso.h" 28 #include "SystemMatrix.h" 29 #include "Solver.h" 30 #include "escript/blocktimer.h" 31 32 /***********************************************************************************/ 33 34 /* free space */ 35 36 void Paso_Solver_free(Paso_SystemMatrix* A) { 37 Paso_Preconditioner_free(A->solver); 38 A->solver=NULL; 39 } 40 /* call the iterative solver: */ 41 42 void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options,Paso_Performance* pp) { 43 double norm2_of_b,tol,tolerance,time_iter,*r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b, 44 norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local,norm_max_of_residual_local,norm_max_of_residual, 45 last_norm_max_of_residual,*scaling, loc_norm; 46 dim_t i,totIter,cntIter,method; 47 bool_t finalizeIteration; 48 err_t errorCode; 49 dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A); 50 dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A); 51 double blocktimer_start = blocktimer_time(); 52 53 tolerance=MAX(options->tolerance,EPSILON); 54 Paso_resetError(); 55 method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric); 56 /* check matrix type */ 57 if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) { 58 Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0."); 59 } 60 if (A->col_block_size != A->row_block_size) { 61 Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal."); 62 } 63 if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) { 64 Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix."); 65 return; 66 } 67 Performance_startMonitor(pp,PERFORMANCE_ALL); 68 if (Paso_noError()) { 69 /* get normalization */ 70 scaling=Paso_SystemMatrix_borrowNormalization(A); 71 if (Paso_noError()) { 72 /* get the norm of the right hand side */ 73 norm2_of_b=0.; 74 norm_max_of_b=0.; 75 #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local) 76 { 77 norm2_of_b_local=0.; 78 norm_max_of_b_local=0.; 79 #pragma omp for private(i) schedule(static) 80 for (i = 0; i < numEqua ; ++i) { 81 norm2_of_b_local += b[i] * b[i]; 82 norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local); 83 } 84 #pragma omp critical 85 { 86 norm2_of_b += norm2_of_b_local; 87 norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b); 88 } 89 } 90 /* TODO: use one call */ 91 #ifdef PASO_MPI 92 { 93 loc_norm = norm2_of_b; 94 MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm); 95 loc_norm = norm_max_of_b; 96 MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm); 97 } 98 #endif 99 norm2_of_b=sqrt(norm2_of_b); 100 101 /* if norm2_of_b==0 we are ready: x=0 */ 102 if (norm2_of_b <=0.) { 103 #pragma omp parallel for private(i) schedule(static) 104 for (i = 0; i < numSol; i++) x[i]=0.; 105 if (options->verbose) printf("right hand side is identical zero.\n"); 106 } else { 107 if (options->verbose) { 108 printf("l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b); 109 printf("l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance); 110 switch (method) { 111 case PASO_BICGSTAB: 112 printf("Iterative method is BiCGStab.\n"); 113 break; 114 case PASO_PCG: 115 printf("Iterative method is PCG.\n"); 116 break; 117 case PASO_PRES20: 118 printf("Iterative method is PRES20.\n"); 119 break; 120 case PASO_GMRES: 121 if (options->restart>0) { 122 printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart); 123 } else { 124 printf("Iterative method is GMRES(%d)\n",options->truncation); 125 } 126 break; 127 } 128 } 129 /* construct the preconditioner */ 130 131 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT); 132 Paso_Solver_setPreconditioner(A,options); 133 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT); 134 if (! Paso_noError()) return; 135 136 time_iter=Paso_timer(); 137 Paso_SystemMatrix_allocBuffer(A); 138 /* get an initial guess by evaluating the preconditioner */ 139 #pragma omp parallel 140 { 141 Paso_Solver_solvePreconditioner(A,x,b); 142 } 143 /* start the iteration process :*/ 144 r=TMPMEMALLOC(numEqua,double); 145 Paso_checkPtr(r); 146 if (Paso_noError()) { 147 totIter = 0; 148 finalizeIteration = FALSE; 149 last_norm2_of_residual=norm2_of_b; 150 last_norm_max_of_residual=norm_max_of_b; 151 while (! finalizeIteration) { 152 cntIter = options->iter_max - totIter; 153 finalizeIteration = TRUE; 154 /* Set initial residual. */ 155 norm2_of_residual = 0; 156 norm_max_of_residual = 0; 157 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local) 158 { 159 #pragma omp for private(i) schedule(static) 160 for (i = 0; i < numEqua; i++) { 161 r[i]=b[i]; 162 } 163 164 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r); 165 166 norm2_of_residual_local = 0; 167 norm_max_of_residual_local = 0; 168 #pragma omp for private(i) schedule(static) 169 for (i = 0; i < numEqua; i++) { 170 norm2_of_residual_local+= r[i] * r[i]; 171 norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local); 172 } 173 #pragma omp critical 174 { 175 norm2_of_residual += norm2_of_residual_local; 176 norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual); 177 } 178 } 179 /* TODO: use one call */ 180 #ifdef PASO_MPI 181 { 182 loc_norm = norm2_of_residual; 183 MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm); 184 loc_norm = norm_max_of_residual; 185 MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm); 186 } 187 #endif 188 norm2_of_residual =sqrt(norm2_of_residual); 189 190 if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual); 191 if (totIter>0 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) { 192 if (options->verbose) printf(" divergence!\n"); 193 Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up."); 194 } else { 195 /* */ 196 if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) { 197 tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b); 198 if (options->verbose) printf(" (new tolerance = %e).\n",tol); 199 last_norm2_of_residual=norm2_of_residual; 200 last_norm_max_of_residual=norm_max_of_residual; 201 /* call the solver */ 202 switch (method) { 203 case PASO_BICGSTAB: 204 errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp); 205 break; 206 case PASO_PCG: 207 errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp); 208 break; 209 case PASO_PRES20: 210 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp); 211 break; 212 case PASO_GMRES: 213 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp); 214 break; 215 } 216 totIter += cntIter; 217 218 /* error handling */ 219 if (errorCode==NO_ERROR) { 220 finalizeIteration = FALSE; 221 } else if (errorCode==SOLVER_MAXITER_REACHED) { 222 Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion."); 223 if (options->verbose) printf("Maximum number of iterations reached.!\n"); 224 } else if (errorCode == SOLVER_INPUT_ERROR ) { 225 Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver."); 226 if (options->verbose) printf("Internal error!\n"); 227 } else if ( errorCode == SOLVER_BREAKDOWN ) { 228 if (cntIter <= 1) { 229 Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver."); 230 if (options->verbose) printf("Uncurable break down!\n"); 231 } else { 232 if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol); 233 finalizeIteration = FALSE; 234 } 235 } else if (errorCode == SOLVER_MEMORY_ERROR) { 236 Paso_setError(MEMORY_ERROR,"memory allocation failed."); 237 if (options->verbose) printf("Memory allocation failed!\n"); 238 } else if (errorCode !=SOLVER_NO_ERROR ) { 239 Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver."); 240 if (options->verbose) printf("Unidentified error!\n"); 241 } 242 } else { 243 if (options->verbose) printf(". convergence! \n"); 244 } 245 } 246 } /* while */ 247 } 248 MEMFREE(r); 249 Paso_SystemMatrix_freeBuffer(A); 250 time_iter=Paso_timer()-time_iter; 251 if (options->verbose) { 252 printf("timing: solver: %.4e sec\n",time_iter); 253 if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter); 254 } 255 } 256 } 257 } 258 Performance_stopMonitor(pp,PERFORMANCE_ALL); 259 blocktimer_increment("Paso_Solver()", blocktimer_start); 260 }

Properties

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