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

Revision 633 - (hide annotations)
Thu Mar 23 05:37:00 2006 UTC (14 years ago) by dhawcroft
Original Path: trunk/paso/src/Solvers/Solver.c
File MIME type: text/plain
File size: 11607 byte(s)

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

Properties

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