00001 /* 00002 GMRES - RESiduum minimization based iterative linear solver 00003 00004 * Written by L. Weimann 00005 * Purpose Iterative solution of large linear systems 00006 * Category ???. - Linear systems 00007 * Keywords large linear system, iterative solver 00008 * Version 1.0 00009 * Revision June 2006 00010 * Latest Change June 2006 00011 * Library NewtonLib 00012 * Code C, Double Precision 00013 * Environment Standard C environment on PC's, 00014 workstations and hosts. 00015 * Copyright (c) Konrad-Zuse-Zentrum fuer 00016 Informationstechnik Berlin (ZIB) 00017 Takustrasse 7, D-14195 Berlin-Dahlem 00018 phone : + 49/30/84185-0 00019 fax : + 49/30/84185-125 00020 * Contact Lutz Weimann 00021 ZIB, Division Scientific Computing, 00022 Department Numerical Analysis and Modelling 00023 phone : + 49/30/84185-185 00024 fax : + 49/30/84185-107 00025 e-mail: weimann@zib.de 00026 00027 * References: 00028 00029 P. Deuflhard: 00030 Newton Methods for Nonlinear Problems. - 00031 Affine Invariance and Adaptive Algorithms. 00032 Series Computational Mathematics 35, Springer (2004) 00033 00034 --------------------------------------------------------------- 00035 00036 * Licence 00037 You may use or modify this code for your own non commercial 00038 purposes for an unlimited time. 00039 In any case you should not deliver this code without a special 00040 permission of ZIB. 00041 In case you intend to use the code commercially, we oblige you 00042 to sign an according licence agreement with ZIB. 00043 00044 * Warranty 00045 This code has been tested up to a certain level. Defects and 00046 weaknesses, which may be included in the code, do not establish 00047 any warranties by ZIB. ZIB does not take over any liabilities 00048 which may follow from acquisition or application of this code. 00049 00050 * Software status 00051 This code is under care of ZIB and belongs to ZIB software class 2. 00052 00053 ------------------------------------------------------------ 00054 00055 * Parameters description 00056 ====================== 00057 00058 The calling interface looks as follows: 00059 00060 extern void gmres(int n, double *y, MATVEC *matvec, 00061 PRECON *preconr, PRECON *preconl, double *b, 00062 struct ITLIN_OPT *opt, struct ITLIN_INFO *info) 00063 00064 The structures used within the parameter list are defined 00065 as follows: 00066 --- 00067 struct ITLIN_OPT 00068 { 00069 double tol, rho; 00070 int i_max, maxiter; 00071 TERM_CHECK termcheck; /* GMRES only * / 00072 CONV_CHECK convcheck; /* PCG only * / 00073 LOGICAL rescale; /* GBIT only * / 00074 PRINT_LEVEL errorlevel, monitorlevel, datalevel; 00075 FILE *errorfile, *monitorfile, *datafile, 00076 *iterfile, *resfile, *miscfile; 00077 double *scale; 00078 }; 00079 00080 where the applicable types used within this structure are defined by 00081 typedef enum {None=0, Minimum=1, Verbose=2, Debug=3} PRINT_LEVEL; 00082 typedef enum { CheckOnRestart=0, CheckEachIter=1 } TERM_CHECK ; 00083 --- 00084 struct ITLIN_INFO 00085 { 00086 double precision, normdx, residuum; 00087 int iter, rcode, subcode, nomatvec, noprecon, noprecl, noprecr; 00088 }; 00089 --- 00090 00091 A detailed description of the parameters follows: 00092 00093 int n : 00094 The number of equations and unknown variables of the linear system. 00095 00096 double *y : 00097 A pointer to an array of double values of size n. 00098 The array must contain on input an initial guess of the linear system 00099 solution, which is used as the start-vector of the iteration. 00100 On output, the pointed array contains an approximate solution vector y*, 00101 which fits the preconditioned residuum reduction condition 00102 ||B_l*r^i||/||B_l*r^0|| <= opt->tol, 00103 where ||w|| denotes the Euclidian norm of w, and B_l denotes the matrix 00104 of the left preconditioner. 00105 00106 void *matvec(int n, double *y, double *z); 00107 A pointer to the matrix times vector multiplication user routine. 00108 This routine is required - no default routine is supplied. 00109 The parameters are: 00110 int n input Number of vector components. 00111 double *y input Vector of unknowns, of size n . 00112 double *z output Vector which holds the matrix-vector product A*y. 00113 00114 void *preconr(int n, double *z, double *w); 00115 A pointer to the right preconditioner user routine, which computes 00116 w=B_r*z, where B_r should be an approximation of the inverse of the 00117 matrix A. If a null pointer is supplied, then a dummy preconditioner 00118 routine will be used which realizes the preconditioner matrix 00119 B_r=identity. 00120 int n input Number of vector components. 00121 double *z input Preconditioned iterate, of size n . 00122 double *w output Vector which holds the matrix-vector product B_r*z, 00123 i.e. the original iterate. 00124 00125 void *preconl(int n, double *z, double *w); 00126 A pointer to the left preconditioner user routine, which computes 00127 w=B_l*z, where B_l should be an approximation of the inverse of the 00128 matrix A. If a null pointer is supplied, then a dummy preconditioner 00129 routine will be used which realizes the preconditioner matrix 00130 B_l=identity. 00131 int n input Number of vector components. 00132 double *z input Residual vector, of size n . 00133 double *w output Vector which holds the matrix-vector product B_l*z, 00134 i.e. the preconditioned residuum. 00135 00136 double *b : 00137 A pointer to an array of double values of size n. 00138 The pointed array must hold the right hand side of the linear system 00139 to solve. 00140 00141 struct ITLIN_OPT *opt: 00142 A pointer to an options structure. The pointed fields of the structure 00143 contain input options to GMRES. 00144 00145 opt->tol is of type double and must contain the error threshold 00146 which the (left preconditioned) residuum norm reduction quotient must fit. 00147 00148 opt->maxiter is of type int and must contain the maximum number of allowed 00149 iterations. if a zero or negative value is supplied, then the maximum 00150 iteration count will be set to 100. 00151 00152 opt->i_max is of type int and must contain the maximum number of 00153 iterations before a restart occurs. The main portion of used memory 00154 by GMRES depends on i_max, such that n*i_max elements of double 00155 storage will be used. If a nonpositive value is supplied, then i_max 00156 is set to 10. 00157 00158 opt->termcheck is of type TERM_CHECK. 00159 If set to CheckOnRestart, then the residuum norm reduction quotient 00160 will be only checked when a restart occurs - such saving additional 00161 computation effort which is necessary on each intermediate iteration 00162 step to compute the quantity ||B_l*A*B_r*y^i||: One preconr call, 00163 one matvec call, and one preconl call. 00164 If set to CheckEachIter, then additional computations on each iteration 00165 will be done to compute the (left preconditioned) residuum. This roughly 00166 doubles the number of calls to the routines matvec, preconr and preconl. 00167 For additional info, refer to the description of the parameter *y above. 00168 00169 opt->errorlevel is of type PRINT_LEVEL. If it is set to level None, 00170 then no error message will be printed if an error condition occurs. 00171 If it is set to level Minimum or any higher level, then error messages 00172 will be printed, if appropriate. 00173 00174 opt->monitorlevel is of type PRINT_LEVEL. If it is set to level None, 00175 then no monitor output will be printed. 00176 If it is set to level Minimum, a few infomation will be printed. 00177 If set to level Verbose, then some infomation about each iteration 00178 step, fitting into a single line, will be printed. This only applies, 00179 whenopt->termcheck=CheckEachIter is set, otherwise information will 00180 be only printed out when a restart occurs. The higher level Debug 00181 is reserved for future additional information output. 00182 00183 opt->datalevel is of type PRINT_LEVEL. If it is set to level None, 00184 then no data output will be printed. 00185 If it is set to level Minimum, then the values of the initial iteration 00186 vector x and the final vector x will be printed. 00187 If set to level Verbose, then the iteration vector x will be printed for 00188 each step. The higher level Debug is reserved for future additional 00189 information output. 00190 00191 opt->errorfile is of type pointer to FILE, as defined by the <stdio.h> 00192 header file. It must be set either to a NULL pointer, stderr, stdout, 00193 or to another file pointer which has been initialized by a fopen call. 00194 If it is set to NULL, opt->errorfile will be set to stdout. The error 00195 messages will be printed to opt->errorfile. 00196 00197 opt->monitorfile is of type pointer to FILE, as defined by the <stdio.h> 00198 header file. It must be set either to a NULL pointer, stderr, stdout, 00199 or to another file pointer which has been initialized by a fopen call. 00200 If it is set to NULL, opt->monitorfile will be set to stdout. The monitor 00201 output will be printed to opt->monitorfile. 00202 00203 opt->datafile is of type pointer to FILE, as defined by the <stdio.h> 00204 header file. It must be set either to a NULL pointer, stderr, stdout, 00205 or to another file pointer which has been initialized by a fopen call. 00206 If it is set to NULL, a file named "gmres.data" will be opened by a 00207 fopen call and opt->datafile will be set to the filepointer which the 00208 fopen returns. The data output will be printed to opt->datafile. 00209 00210 opt->iterfile is of type pointer to FILE, as defined by the <stdio.h> 00211 header file. It must be set either to a NULL pointer or to file pointer 00212 which has been initialized by a fopen call. The iteration number and 00213 the iteration vector will be written out to the associated file, for 00214 each iteration step. If opt->iterfile is set to NULL, no such 00215 data will be written out. 00216 00217 opt->resfile is of type pointer to FILE, as defined by the <stdio.h> 00218 header file. It must be set either to a NULL pointer or to file pointer 00219 which has been initialized by a fopen call. The iteration number and 00220 the residuum vector will be written out to the associated file, for 00221 each iteration step. If opt->resfile is set to NULL, no such 00222 data will be written out. 00223 00224 opt->miscfile is of type pointer to FILE, as defined by the <stdio.h> 00225 header file. It must be set either to a NULL pointer or to file pointer 00226 which has been initialized by a fopen call. The iteration number, an 00227 identification number of the calling code (0 for GMRES), the norm of the 00228 (left preconditioned) residuum, followed by three zero floating point 00229 values as placeholders, will be written out, for each iteration step, 00230 if opt->termcheck=CheckEachIter is set. Otherwise, if 00231 opt->termcheck=CheckOnRestart is set, data will be only written out 00232 when a restart occurs. 00233 If opt->miscfile is set to NULL, no such data will be written out. 00234 00235 Note: The output to the files opt->iterfile, opt->resfile and 00236 opt->miscfile is written as a single for each iteration step. 00237 Such, the data in these files are suitable as input to the 00238 graphics utility GNUPLOT. 00239 00240 struct ITLIN_INFO *info: 00241 A pointer to an info structure. The pointed fields of the structure 00242 are set output info of GMRES. 00243 00244 info->precision is of type double and is set to the achieved residual 00245 reduction ||r^i||/||r^0||. 00246 00247 info->iter is set to number of iteration steps done. 00248 00249 info->nomatvec is set to the number of done calls to the matrix times 00250 vector multiplication user routine matvec. 00251 00252 info->noprecr is set to the number of done calls to the right 00253 preconditioner user routine preconr or the dummy preconditioner routine, 00254 if the user didn't supply a right preconditioner routine. 00255 00256 info->noprecl is set to the number of done calls to the left 00257 preconditioner user routine preconl or the dummy preconditioner routine, 00258 if the user didn't supply a left preconditioner routine. 00259 00260 info->rcode is set to the return-code of GMRES. A return-code 0 00261 means that GMRES has terminated sucessfully. For the meaning of a 00262 nonzero return-code, see the error messages list below. 00263 00264 00265 The following error conditions may occur: (returned via info->rcode) 00266 -------------------------------------------------------------------- 00267 00268 -999 routine zibnum_fwalloc failed to allocate double memory via malloc. 00269 -997 routine zibnum_pfwalloc failed to allocate double pointer memory 00270 via malloc. 00271 -995 Internal i/o control block could not be allocated via malloc. 00272 -994 Internally used data structure could not be allocated via malloc. 00273 -989 Default data-output file could not be opened via fopen call. 00274 -99 NULL pointer supplied for matvec - the matrix times vector routine 00275 must be defined! 00276 2 Maximum number of iterations (as set by opt->maxiter) exceeded. 00277 20 Nonpositive input for dimensional parameter n. 00278 21 Nonpositive value for opt->tol supplied. 00279 00280 Summary of changes: 00281 ------------------- 00282 00283 Version Date Changes 00284 0.1 2006/06/13 Initial Prerelease. 00285 00286 */ 00287 #include <stdlib.h> 00288 #include <stdio.h> 00289 #include <math.h> 00290 #include "itlin.h" 00291 00292 double gmres_norm2(int n, double *v); 00293 int gmres_qrfact(int n, int lda, double *a, double *q, int ijob); 00294 void gmres_qrsolv(int n, int lda, double *a, double *q, double *b); 00295 00296 #define MAXITER_DEFAULT 100 00297 00298 extern struct ITLIN_IO *itlin_ioctl; 00299 00300 extern void gmres(int n, double *y, MATVEC *matvec, 00301 PRECON *preconr, PRECON *preconl, double *b, 00302 struct ITLIN_OPT *opt, struct ITLIN_INFO *info) 00303 { 00304 int i, j, l, m, m1, m11, l1, iter=0, fail,ijob, 00305 nomatvec=0, nopreconl=0, 00306 nopreconr=0, i_max, max_iter = opt->maxiter; 00307 double **V, *r, *y0, *w, *vi, *vip1, *vl, *Hi, *h, *z, *q; 00308 double s, normvip1, beta, beta0, beta1, eta, tol=opt->tol; 00309 LOGICAL stop_iter, io_allocated=False; 00310 TERM_CHECK termcheck = opt->termcheck; 00311 struct ITLIN_DATA *data=malloc(sizeof(struct ITLIN_DATA)); 00312 00313 if (!itlin_ioctl) itlin_ioctl=malloc(sizeof(struct ITLIN_IO)); 00314 if (!itlin_ioctl) 00315 { fprintf(stderr,"\n could not allocate output controlblock\n"); 00316 RCODE=-995; return; } 00317 else 00318 io_allocated = True; 00319 if (!data) 00320 { fprintf(stderr,"\n could not allocate struct data\n"); 00321 RCODE=-994; return; }; 00322 data->codeid = GMRES; 00323 data->normdx = 0.0; 00324 data->tau = 0.0; 00325 data->t = 0.0; 00326 data->mode = Initial; 00327 ERRORLEVEL = opt->errorlevel; 00328 MONITORLEVEL = opt->monitorlevel; 00329 DATALEVEL = opt->datalevel; 00330 ERROR = opt->errorfile; 00331 MONITOR = opt->monitorfile; 00332 DATA = opt->datafile; 00333 FITER = opt->iterfile; 00334 FRES = opt->resfile; 00335 FMISC = opt->miscfile; 00336 if ( !ERROR && ERRORLEVEL>0 ) ERROR = stdout; 00337 if ( !MONITOR && MONITORLEVEL>0 ) MONITOR = stdout; 00338 if ( !DATA && DATALEVEL>0 ) 00339 { DATA=fopen("gmres.data","w"); 00340 if (!DATA && ERRORLEVEL>0) 00341 { fprintf(ERROR,"\n fopen of file gmres.data failed\n"); 00342 RCODE=-989; return; 00343 }; 00344 }; 00345 opt->errorfile = ERROR; 00346 opt->monitorfile = MONITOR; 00347 opt->datafile = DATA; 00348 if ( MONITORLEVEL > 0 ) fprintf(MONITOR,"\n GMRES - Version 0.1\n"); 00349 if ( max_iter <= 0 ) max_iter = MAXITER_DEFAULT; 00350 RCODE = itlin_parcheck_and_print(n,matvec,opt,0); 00351 if ( RCODE !=0 ) 00352 { if (io_allocated) {free(itlin_ioctl); itlin_ioctl=NULL;}; 00353 if (data) free(data); 00354 return; 00355 }; 00356 i_max = opt->i_max; 00357 if (!preconr) preconr = &itlin_noprecon; 00358 if (!preconl) preconl = &itlin_noprecon; 00359 RCODE = zibnum_pfwalloc(i_max+2,&V,"V"); if ( RCODE !=0 ) return; 00360 RCODE = zibnum_fwalloc(i_max+1,&h,"h"); if ( RCODE !=0 ) return; 00361 RCODE = zibnum_fwalloc(i_max+1,&z,"z"); if ( RCODE !=0 ) return; 00362 RCODE = zibnum_fwalloc((i_max+1)*(i_max+1),&Hi,"Hi"); 00363 if ( RCODE !=0 ) return; 00364 RCODE = zibnum_fwalloc(n,&r,"r"); if ( RCODE !=0 ) return; 00365 RCODE = zibnum_fwalloc(n,&y0,"y0"); if ( RCODE !=0 ) return; 00366 RCODE = zibnum_fwalloc(n,&w,"w"); if ( RCODE !=0 ) return; 00367 RCODE = zibnum_fwalloc(2*(i_max+1),&q,"q"); if ( RCODE !=0 ) return; 00368 for (i=0;i<=i_max+1;i++) 00369 { RCODE = zibnum_fwalloc(n,&V[i],"V[i]"); 00370 if ( RCODE !=0 ) return; 00371 }; 00372 /* Initialization */ 00373 if ( MONITORLEVEL>1 ) 00374 fprintf(MONITOR,"\n\n iter i norm(res) eta\n\n"); 00375 data->res = r; 00376 restart: 00377 if ( iter > 0 && MONITORLEVEL > 1 && termcheck==CheckEachIter ) 00378 fprintf(MONITOR," > RESTART\n"); 00379 for (j=0;j<n;j++) y0[j]=y[j]; 00380 preconr(n,y,V[0]); nopreconr++; 00381 matvec(n,V[0],w); nomatvec++; 00382 preconl(n,w,r); nopreconl++; 00383 for (j=0;j<n;j++) r[j] = b[j]-r[j]; 00384 beta = gmres_norm2(n,r); 00385 if ( iter==0 ) 00386 { beta0 = beta; eta = ( beta0==0.0 ? 0.0 : 1.0 ); } 00387 else 00388 eta = beta/beta0; 00389 beta1 = beta; 00390 vi = V[0]; 00391 for (j=0;j<n;j++) vi[j] = r[j]/beta; 00392 00393 for (i=1; i<=i_max && eta>tol; i++) 00394 { 00395 if ( i==1 || termcheck==CheckEachIter ) 00396 { 00397 data->residuum = beta1; itlin_dataout(iter,n,y,data); 00398 data->mode = Intermediate; 00399 if ( MONITORLEVEL>1 ) 00400 fprintf(MONITOR," %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta); 00401 }; 00402 00403 /* I. Orthogonalization */ 00404 vi = V[i-1]; 00405 preconr(n,vi,V[i]); nopreconr++; 00406 matvec(n,V[i],w); nomatvec++; 00407 preconl(n,w,r); nopreconl++; 00408 for (l=0;l<i;l++) 00409 { vl = V[l]; s=0.0; 00410 for (j=0;j<n;j++) s += vl[j]*r[j]; h[l] = s; 00411 }; 00412 vip1 = V[i]; 00413 for (j=0;j<n;j++) 00414 { s = 0.0; for (l=0;l<i;l++) s += V[l][j]*h[l]; 00415 vip1[j] = r[j]-s; 00416 }; 00417 00418 /* II. Normalization */ 00419 normvip1 = gmres_norm2(n,vip1); 00420 for (j=0;j<n;j++) vip1[j] /= normvip1; 00421 00422 /* III. Update */ 00423 if ( i>1 ) 00424 { for (l=0;l<i;l++) Hi[(i_max+1)*(i-1)+l] = h[l]; 00425 for (l=0;l<i-1;l++) Hi[(i_max+1)*l+i] = 0.0; 00426 Hi[(i_max+1)*(i-1)+i] = normvip1; 00427 } 00428 else 00429 { 00430 Hi[0] = h[0]; Hi[1] = normvip1; 00431 }; 00432 00433 /* IV. Least squares problem for z*/ 00434 ijob = ( i==1 ? 1 : 2); 00435 fail = gmres_qrfact(i,i_max+1,Hi,q,ijob); 00436 if ( fail != 0 ) 00437 { RCODE = 1; fprintf(ERROR,"\n gmres_qrfact failed with code=%i\n",fail); 00438 goto errorexit; 00439 }; 00440 if ( termcheck==CheckEachIter || i==i_max ) 00441 { 00442 for (l=0;l<i+1;l++) z[l]=0.0; z[0]=beta; 00443 gmres_qrsolv(i,i_max+1,Hi,q,z); 00444 00445 /* V. Approximate solution*/ 00446 for (j=0;j<n;j++) 00447 { s=0.0; 00448 for (l=0;l<i;l++) s+= V[l][j]*z[l]; 00449 y[j] = y0[j]+s; 00450 }; 00451 if ( termcheck==CheckEachIter ) 00452 { 00453 preconr(n,y,V[i+1]); nopreconr++; 00454 matvec(n,V[i+1],w); nomatvec++; 00455 preconl(n,w,r); nopreconl++; 00456 for (j=0;j<n;j++) r[j] = b[j]-r[j]; 00457 beta1 = gmres_norm2(n,r); 00458 eta = beta1/beta0; 00459 }; 00460 }; 00461 iter++; 00462 }; 00463 00464 if ( termcheck==CheckEachIter && MONITORLEVEL>1 ) 00465 fprintf(MONITOR," %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta); 00466 if ( eta > tol && iter < max_iter ) goto restart; 00467 else RCODE = ( eta <= tol ? 0 : 2 ); 00468 if ( MONITORLEVEL>1 && termcheck==CheckOnRestart ) 00469 fprintf(MONITOR," %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta); 00470 data->mode = ( RCODE == 0 ? Solution : Final ); 00471 data->residuum = beta1; itlin_dataout(iter,n,y,data); 00472 00473 errorexit: 00474 if ( ERRORLEVEL > 0 && RCODE != 0 ) 00475 { 00476 switch ( RCODE ) 00477 { 00478 case 2: 00479 fprintf(ERROR,"\n Error - Maximum allowed number of iterations exceeded\n"); 00480 break; 00481 default : 00482 fprintf(ERROR,"\n Error, code=%i\n",RCODE); 00483 }; 00484 }; 00485 for (i=0;i<=i_max+1;i++) free(V[i]); 00486 free(V); free(h); free(z); free(Hi); free(r); free(y0); 00487 00488 info->precision = eta; 00489 info->iter = iter; 00490 info->nomatvec = nomatvec; 00491 info->noprecl = nopreconl; 00492 info->noprecr = nopreconr; 00493 } 00494 00495 double gmres_norm2(int n, double *v) 00496 { int i; 00497 double rval = 0.0; 00498 for (i=0;i<n;i++) rval += v[i]*v[i]; 00499 return sqrt( rval ); 00500 } 00501 00502 int gmres_qrfact(int n, int lda, double *a, double *q, int ijob) 00503 /* 00504 Translation of Fortran routine DHEQR from file DGMRES of 00505 library=slatec(slap) package. 00506 */ 00507 { int i, j, j1, k,info=0; 00508 double t, t1, t2, c, s; 00509 if ( ijob <= 1 ) 00510 { for (k=0;k<n;k++) 00511 { for (j=0;j<k-1;j++) 00512 { i=2*j; c=q[i]; s=q[i+1]; 00513 j1=j+k*lda; t1=a[j1]; t2=a[j1+1]; 00514 a[j1]=c*t1-s*t2; a[j1+1]=s*t1+c*t2; 00515 }; 00516 i=2*k; j1=k+lda*k; t1=a[j1]; t2=a[j1+1]; 00517 if ( t2==0.0 ) { c=1.0; s=0.0; } 00518 else if ( fabs(t2)>=fabs(t1) ) 00519 { t=t1/t2; s=-1.0/sqrt(1.0+t*t); c=-s*t; } 00520 else 00521 { t=t2/t1; c=1.0/sqrt(1.0+t*t); s=-c*t; }; 00522 q[i]=c; q[i+1]=s; a[j1]=c*t1-s*t2; 00523 if( a[j1]==0.0 ) info = k; 00524 }; 00525 } 00526 else 00527 { for (k=0;k<n-1;k++) 00528 { 00529 i=2*k; j1=k+lda*(n-1); 00530 t1=a[j1]; t2=a[j1+1]; 00531 c=q[i]; s=q[i+1]; 00532 a[j1]=c*t1-s*t2; a[j1+1]=s*t1 + c*t2; 00533 }; 00534 info=0; 00535 j1 = n-1+lda*(n-1); t1 = a[j1]; t2 = a[j1+1]; 00536 if ( t2==0.0 ) { c=1.0; s=0.0; } 00537 else if ( fabs(t2)>=fabs(t1) ) 00538 { t=t1/t2; s=-1.0/sqrt(1.0+t*t); c=-s*t; } 00539 else 00540 { t=t2/t1; c=1.0/sqrt(1.0+t*t); s=-c*t; }; 00541 i=2*n-2; q[i]=c; q[i+1]=s; a[j1]=c*t1-s*t2; 00542 if( a[j1]==0.0 ) info = n-1; 00543 }; 00544 return info; 00545 } 00546 00547 void gmres_qrsolv(int n, int lda, double *a, double *q, double *b) 00548 /* 00549 Translation of Fortran routine DHELS from file DGMRES of 00550 library=slatec(slap) package. 00551 */ 00552 { int i, k, kb; 00553 double t, t1, t2, c, s; 00554 for (k=0;k<n;k++) 00555 { i=2*k; c=q[i]; s=q[i+1]; 00556 t1=b[k]; t2=b[k+1]; 00557 b[k]=c*t1-s*t2; b[k+1]=s*t1+c*t2; 00558 }; 00559 for (kb=0;kb<n;kb++) 00560 { 00561 k=n-1-kb; b[k] /= a[k+lda*k]; 00562 t = -b[k]; 00563 for (i=0;i<k;i++) b[i] += t*a[i+k*lda]; 00564 }; 00565 return; 00566 }