VERB_code_2.3
gmres.c
1 /*
2  GMRES - RESiduum minimization based iterative linear solver
3 
4  * Written by L. Weimann
5  * Purpose Iterative solution of large linear systems
6  * Category ???. - Linear systems
7  * Keywords large linear system, iterative solver
8  * Version 1.0
9  * Revision June 2006
10  * Latest Change June 2006
11  * Library NewtonLib
12  * Code C, Double Precision
13  * Environment Standard C environment on PC's,
14  workstations and hosts.
15  * Copyright (c) Konrad-Zuse-Zentrum fuer
16  Informationstechnik Berlin (ZIB)
17  Takustrasse 7, D-14195 Berlin-Dahlem
18  phone : + 49/30/84185-0
19  fax : + 49/30/84185-125
20  * Contact Lutz Weimann
21  ZIB, Division Scientific Computing,
22  Department Numerical Analysis and Modelling
23  phone : + 49/30/84185-185
24  fax : + 49/30/84185-107
25  e-mail: weimann@zib.de
26 
27 * References:
28 
29  P. Deuflhard:
30  Newton Methods for Nonlinear Problems. -
31  Affine Invariance and Adaptive Algorithms.
32  Series Computational Mathematics 35, Springer (2004)
33 
34  ---------------------------------------------------------------
35 
36  * Licence
37  You may use or modify this code for your own non commercial
38  purposes for an unlimited time.
39  In any case you should not deliver this code without a special
40  permission of ZIB.
41  In case you intend to use the code commercially, we oblige you
42  to sign an according licence agreement with ZIB.
43 
44  * Warranty
45  This code has been tested up to a certain level. Defects and
46  weaknesses, which may be included in the code, do not establish
47  any warranties by ZIB. ZIB does not take over any liabilities
48  which may follow from acquisition or application of this code.
49 
50  * Software status
51  This code is under care of ZIB and belongs to ZIB software class 2.
52 
53  ------------------------------------------------------------
54 
55  * Parameters description
56  ======================
57 
58  The calling interface looks as follows:
59 
60  extern void gmres(int n, double *y, MATVEC *matvec,
61  PRECON *preconr, PRECON *preconl, double *b,
62  struct ITLIN_OPT *opt, struct ITLIN_INFO *info)
63 
64  The structures used within the parameter list are defined
65  as follows:
66  ---
67  struct ITLIN_OPT
68  {
69  double tol, rho;
70  int i_max, maxiter;
71  TERM_CHECK termcheck; /* GMRES only * /
72  CONV_CHECK convcheck; /* PCG only * /
73  LOGICAL rescale; /* GBIT only * /
74  PRINT_LEVEL errorlevel, monitorlevel, datalevel;
75  FILE *errorfile, *monitorfile, *datafile,
76  *iterfile, *resfile, *miscfile;
77  double *scale;
78  };
79 
80  where the applicable types used within this structure are defined by
81  typedef enum {None=0, Minimum=1, Verbose=2, Debug=3} PRINT_LEVEL;
82  typedef enum { CheckOnRestart=0, CheckEachIter=1 } TERM_CHECK ;
83  ---
84  struct ITLIN_INFO
85  {
86  double precision, normdx, residuum;
87  int iter, rcode, subcode, nomatvec, noprecon, noprecl, noprecr;
88  };
89  ---
90 
91  A detailed description of the parameters follows:
92 
93  int n :
94  The number of equations and unknown variables of the linear system.
95 
96  double *y :
97  A pointer to an array of double values of size n.
98  The array must contain on input an initial guess of the linear system
99  solution, which is used as the start-vector of the iteration.
100  On output, the pointed array contains an approximate solution vector y*,
101  which fits the preconditioned residuum reduction condition
102  ||B_l*r^i||/||B_l*r^0|| <= opt->tol,
103  where ||w|| denotes the Euclidian norm of w, and B_l denotes the matrix
104  of the left preconditioner.
105 
106  void *matvec(int n, double *y, double *z);
107  A pointer to the matrix times vector multiplication user routine.
108  This routine is required - no default routine is supplied.
109  The parameters are:
110  int n input Number of vector components.
111  double *y input Vector of unknowns, of size n .
112  double *z output Vector which holds the matrix-vector product A*y.
113 
114  void *preconr(int n, double *z, double *w);
115  A pointer to the right preconditioner user routine, which computes
116  w=B_r*z, where B_r should be an approximation of the inverse of the
117  matrix A. If a null pointer is supplied, then a dummy preconditioner
118  routine will be used which realizes the preconditioner matrix
119  B_r=identity.
120  int n input Number of vector components.
121  double *z input Preconditioned iterate, of size n .
122  double *w output Vector which holds the matrix-vector product B_r*z,
123  i.e. the original iterate.
124 
125  void *preconl(int n, double *z, double *w);
126  A pointer to the left preconditioner user routine, which computes
127  w=B_l*z, where B_l should be an approximation of the inverse of the
128  matrix A. If a null pointer is supplied, then a dummy preconditioner
129  routine will be used which realizes the preconditioner matrix
130  B_l=identity.
131  int n input Number of vector components.
132  double *z input Residual vector, of size n .
133  double *w output Vector which holds the matrix-vector product B_l*z,
134  i.e. the preconditioned residuum.
135 
136  double *b :
137  A pointer to an array of double values of size n.
138  The pointed array must hold the right hand side of the linear system
139  to solve.
140 
141  struct ITLIN_OPT *opt:
142  A pointer to an options structure. The pointed fields of the structure
143  contain input options to GMRES.
144 
145  opt->tol is of type double and must contain the error threshold
146  which the (left preconditioned) residuum norm reduction quotient must fit.
147 
148  opt->maxiter is of type int and must contain the maximum number of allowed
149  iterations. if a zero or negative value is supplied, then the maximum
150  iteration count will be set to 100.
151 
152  opt->i_max is of type int and must contain the maximum number of
153  iterations before a restart occurs. The main portion of used memory
154  by GMRES depends on i_max, such that n*i_max elements of double
155  storage will be used. If a nonpositive value is supplied, then i_max
156  is set to 10.
157 
158  opt->termcheck is of type TERM_CHECK.
159  If set to CheckOnRestart, then the residuum norm reduction quotient
160  will be only checked when a restart occurs - such saving additional
161  computation effort which is necessary on each intermediate iteration
162  step to compute the quantity ||B_l*A*B_r*y^i||: One preconr call,
163  one matvec call, and one preconl call.
164  If set to CheckEachIter, then additional computations on each iteration
165  will be done to compute the (left preconditioned) residuum. This roughly
166  doubles the number of calls to the routines matvec, preconr and preconl.
167  For additional info, refer to the description of the parameter *y above.
168 
169  opt->errorlevel is of type PRINT_LEVEL. If it is set to level None,
170  then no error message will be printed if an error condition occurs.
171  If it is set to level Minimum or any higher level, then error messages
172  will be printed, if appropriate.
173 
174  opt->monitorlevel is of type PRINT_LEVEL. If it is set to level None,
175  then no monitor output will be printed.
176  If it is set to level Minimum, a few infomation will be printed.
177  If set to level Verbose, then some infomation about each iteration
178  step, fitting into a single line, will be printed. This only applies,
179  whenopt->termcheck=CheckEachIter is set, otherwise information will
180  be only printed out when a restart occurs. The higher level Debug
181  is reserved for future additional information output.
182 
183  opt->datalevel is of type PRINT_LEVEL. If it is set to level None,
184  then no data output will be printed.
185  If it is set to level Minimum, then the values of the initial iteration
186  vector x and the final vector x will be printed.
187  If set to level Verbose, then the iteration vector x will be printed for
188  each step. The higher level Debug is reserved for future additional
189  information output.
190 
191  opt->errorfile is of type pointer to FILE, as defined by the <stdio.h>
192  header file. It must be set either to a NULL pointer, stderr, stdout,
193  or to another file pointer which has been initialized by a fopen call.
194  If it is set to NULL, opt->errorfile will be set to stdout. The error
195  messages will be printed to opt->errorfile.
196 
197  opt->monitorfile is of type pointer to FILE, as defined by the <stdio.h>
198  header file. It must be set either to a NULL pointer, stderr, stdout,
199  or to another file pointer which has been initialized by a fopen call.
200  If it is set to NULL, opt->monitorfile will be set to stdout. The monitor
201  output will be printed to opt->monitorfile.
202 
203  opt->datafile is of type pointer to FILE, as defined by the <stdio.h>
204  header file. It must be set either to a NULL pointer, stderr, stdout,
205  or to another file pointer which has been initialized by a fopen call.
206  If it is set to NULL, a file named "gmres.data" will be opened by a
207  fopen call and opt->datafile will be set to the filepointer which the
208  fopen returns. The data output will be printed to opt->datafile.
209 
210  opt->iterfile is of type pointer to FILE, as defined by the <stdio.h>
211  header file. It must be set either to a NULL pointer or to file pointer
212  which has been initialized by a fopen call. The iteration number and
213  the iteration vector will be written out to the associated file, for
214  each iteration step. If opt->iterfile is set to NULL, no such
215  data will be written out.
216 
217  opt->resfile is of type pointer to FILE, as defined by the <stdio.h>
218  header file. It must be set either to a NULL pointer or to file pointer
219  which has been initialized by a fopen call. The iteration number and
220  the residuum vector will be written out to the associated file, for
221  each iteration step. If opt->resfile is set to NULL, no such
222  data will be written out.
223 
224  opt->miscfile is of type pointer to FILE, as defined by the <stdio.h>
225  header file. It must be set either to a NULL pointer or to file pointer
226  which has been initialized by a fopen call. The iteration number, an
227  identification number of the calling code (0 for GMRES), the norm of the
228  (left preconditioned) residuum, followed by three zero floating point
229  values as placeholders, will be written out, for each iteration step,
230  if opt->termcheck=CheckEachIter is set. Otherwise, if
231  opt->termcheck=CheckOnRestart is set, data will be only written out
232  when a restart occurs.
233  If opt->miscfile is set to NULL, no such data will be written out.
234 
235  Note: The output to the files opt->iterfile, opt->resfile and
236  opt->miscfile is written as a single for each iteration step.
237  Such, the data in these files are suitable as input to the
238  graphics utility GNUPLOT.
239 
240  struct ITLIN_INFO *info:
241  A pointer to an info structure. The pointed fields of the structure
242  are set output info of GMRES.
243 
244  info->precision is of type double and is set to the achieved residual
245  reduction ||r^i||/||r^0||.
246 
247  info->iter is set to number of iteration steps done.
248 
249  info->nomatvec is set to the number of done calls to the matrix times
250  vector multiplication user routine matvec.
251 
252  info->noprecr is set to the number of done calls to the right
253  preconditioner user routine preconr or the dummy preconditioner routine,
254  if the user didn't supply a right preconditioner routine.
255 
256  info->noprecl is set to the number of done calls to the left
257  preconditioner user routine preconl or the dummy preconditioner routine,
258  if the user didn't supply a left preconditioner routine.
259 
260  info->rcode is set to the return-code of GMRES. A return-code 0
261  means that GMRES has terminated sucessfully. For the meaning of a
262  nonzero return-code, see the error messages list below.
263 
264 
265  The following error conditions may occur: (returned via info->rcode)
266  --------------------------------------------------------------------
267 
268  -999 routine zibnum_fwalloc failed to allocate double memory via malloc.
269  -997 routine zibnum_pfwalloc failed to allocate double pointer memory
270  via malloc.
271  -995 Internal i/o control block could not be allocated via malloc.
272  -994 Internally used data structure could not be allocated via malloc.
273  -989 Default data-output file could not be opened via fopen call.
274  -99 NULL pointer supplied for matvec - the matrix times vector routine
275  must be defined!
276  2 Maximum number of iterations (as set by opt->maxiter) exceeded.
277  20 Nonpositive input for dimensional parameter n.
278  21 Nonpositive value for opt->tol supplied.
279 
280  Summary of changes:
281  -------------------
282 
283  Version Date Changes
284  0.1 2006/06/13 Initial Prerelease.
285 
286 */
287 #include <stdlib.h>
288 #include <stdio.h>
289 #include <math.h>
290 #include "itlin.h"
291 
292 double gmres_norm2(int n, double *v);
293 int gmres_qrfact(int n, int lda, double *a, double *q, int ijob);
294 void gmres_qrsolv(int n, int lda, double *a, double *q, double *b);
295 
296 #define MAXITER_DEFAULT 100
297 
298 extern struct ITLIN_IO *itlin_ioctl;
299 
300 extern void gmres(int n, double *y, MATVEC *matvec,
301  PRECON *preconr, PRECON *preconl, double *b,
302  struct ITLIN_OPT *opt, struct ITLIN_INFO *info)
303 {
304  int i, j, l, m, m1, m11, l1, iter=0, fail,ijob,
305  nomatvec=0, nopreconl=0,
306  nopreconr=0, i_max, max_iter = opt->maxiter;
307  double **V, *r, *y0, *w, *vi, *vip1, *vl, *Hi, *h, *z, *q;
308  double s, normvip1, beta, beta0, beta1, eta, tol=opt->tol;
309  LOGICAL stop_iter, io_allocated=False;
310  TERM_CHECK termcheck = opt->termcheck;
311  struct ITLIN_DATA *data=malloc(sizeof(struct ITLIN_DATA));
312 
313  if (!itlin_ioctl) itlin_ioctl=malloc(sizeof(struct ITLIN_IO));
314  if (!itlin_ioctl)
315  { fprintf(stderr,"\n could not allocate output controlblock\n");
316  RCODE=-995; return; }
317  else
318  io_allocated = True;
319  if (!data)
320  { fprintf(stderr,"\n could not allocate struct data\n");
321  RCODE=-994; return; };
322  data->codeid = GMRES;
323  data->normdx = 0.0;
324  data->tau = 0.0;
325  data->t = 0.0;
326  data->mode = Initial;
327  ERRORLEVEL = opt->errorlevel;
328  MONITORLEVEL = opt->monitorlevel;
329  DATALEVEL = opt->datalevel;
330  ERROR = opt->errorfile;
331  MONITOR = opt->monitorfile;
332  DATA = opt->datafile;
333  FITER = opt->iterfile;
334  FRES = opt->resfile;
335  FMISC = opt->miscfile;
336  if ( !ERROR && ERRORLEVEL>0 ) ERROR = stdout;
337  if ( !MONITOR && MONITORLEVEL>0 ) MONITOR = stdout;
338  if ( !DATA && DATALEVEL>0 )
339  { DATA=fopen("gmres.data","w");
340  if (!DATA && ERRORLEVEL>0)
341  { fprintf(ERROR,"\n fopen of file gmres.data failed\n");
342  RCODE=-989; return;
343  };
344  };
345  opt->errorfile = ERROR;
346  opt->monitorfile = MONITOR;
347  opt->datafile = DATA;
348  if ( MONITORLEVEL > 0 ) fprintf(MONITOR,"\n GMRES - Version 0.1\n");
349  if ( max_iter <= 0 ) max_iter = MAXITER_DEFAULT;
350  RCODE = itlin_parcheck_and_print(n,matvec,opt,0);
351  if ( RCODE !=0 )
352  {
353  if (io_allocated) {
354  free(itlin_ioctl); itlin_ioctl=NULL;
355  };
356  if (data) free(data);
357  return;
358  };
359  i_max = opt->i_max;
360  if (!preconr) preconr = &itlin_noprecon;
361  if (!preconl) preconl = &itlin_noprecon;
362  RCODE = zibnum_pfwalloc(i_max+2,&V,"V"); if ( RCODE !=0 ) return;
363  RCODE = zibnum_fwalloc(i_max+1,&h,"h"); if ( RCODE !=0 ) return;
364  RCODE = zibnum_fwalloc(i_max+1,&z,"z"); if ( RCODE !=0 ) return;
365  RCODE = zibnum_fwalloc((i_max+1)*(i_max+1),&Hi,"Hi");
366  if ( RCODE !=0 ) return;
367  RCODE = zibnum_fwalloc(n,&r,"r"); if ( RCODE !=0 ) return;
368  RCODE = zibnum_fwalloc(n,&y0,"y0"); if ( RCODE !=0 ) return;
369  RCODE = zibnum_fwalloc(n,&w,"w"); if ( RCODE !=0 ) return;
370  RCODE = zibnum_fwalloc(2*(i_max+1),&q,"q"); if ( RCODE !=0 ) return;
371  for (i=0;i<=i_max+1;i++)
372  { RCODE = zibnum_fwalloc(n,&V[i],"V[i]");
373  if ( RCODE !=0 ) return;
374  };
375  /* Initialization */
376  if ( MONITORLEVEL>1 )
377  fprintf(MONITOR,"\n\n iter i norm(res) eta\n\n");
378  data->res = r;
379  restart:
380  if ( iter > 0 && MONITORLEVEL > 1 && termcheck==CheckEachIter )
381  fprintf(MONITOR," > RESTART\n");
382  for (j=0;j<n;j++) y0[j]=y[j];
383  preconr(n,y,V[0]); nopreconr++;
384  matvec(n,V[0],w); nomatvec++;
385  preconl(n,w,r); nopreconl++;
386  for (j=0;j<n;j++) r[j] = b[j]-r[j];
387  beta = gmres_norm2(n,r);
388  if ( iter==0 )
389  { beta0 = beta; eta = ( beta0==0.0 ? 0.0 : 1.0 ); }
390  else
391  eta = beta/beta0;
392  beta1 = beta;
393  vi = V[0];
394  for (j=0;j<n;j++) vi[j] = r[j]/beta;
395 
396  for (i=1; i<=i_max && eta>tol; i++)
397  {
398  if ( i==1 || termcheck==CheckEachIter )
399  {
400  data->residuum = beta1; itlin_dataout(iter,n,y,data);
401  data->mode = Intermediate;
402  if ( MONITORLEVEL>1 )
403  fprintf(MONITOR," %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta);
404  };
405 
406  /* I. Orthogonalization */
407  vi = V[i-1];
408  preconr(n,vi,V[i]); nopreconr++;
409  matvec(n,V[i],w); nomatvec++;
410  preconl(n,w,r); nopreconl++;
411  for (l=0;l<i;l++)
412  { vl = V[l]; s=0.0;
413  for (j=0;j<n;j++) s += vl[j]*r[j]; h[l] = s;
414  };
415  vip1 = V[i];
416  for (j=0;j<n;j++)
417  { s = 0.0; for (l=0;l<i;l++) s += V[l][j]*h[l];
418  vip1[j] = r[j]-s;
419  };
420 
421  /* II. Normalization */
422  normvip1 = gmres_norm2(n,vip1);
423  for (j=0;j<n;j++) vip1[j] /= normvip1;
424 
425  /* III. Update */
426  if ( i>1 )
427  { for (l=0;l<i;l++) Hi[(i_max+1)*(i-1)+l] = h[l];
428  for (l=0;l<i-1;l++) Hi[(i_max+1)*l+i] = 0.0;
429  Hi[(i_max+1)*(i-1)+i] = normvip1;
430  }
431  else
432  {
433  Hi[0] = h[0]; Hi[1] = normvip1;
434  };
435 
436  /* IV. Least squares problem for z*/
437  ijob = ( i==1 ? 1 : 2);
438  fail = gmres_qrfact(i,i_max+1,Hi,q,ijob);
439  if ( fail != 0 )
440  { RCODE = 1; fprintf(ERROR,"\n gmres_qrfact failed with code=%i\n",fail);
441  goto errorexit;
442  };
443  if ( termcheck==CheckEachIter || i==i_max )
444  {
445  for (l=0;l<i+1;l++) z[l]=0.0; z[0]=beta;
446  gmres_qrsolv(i,i_max+1,Hi,q,z);
447 
448  /* V. Approximate solution*/
449  for (j=0;j<n;j++)
450  { s=0.0;
451  for (l=0;l<i;l++) s+= V[l][j]*z[l];
452  y[j] = y0[j]+s;
453  };
454  if ( termcheck==CheckEachIter )
455  {
456  preconr(n,y,V[i+1]); nopreconr++;
457  matvec(n,V[i+1],w); nomatvec++;
458  preconl(n,w,r); nopreconl++;
459  for (j=0;j<n;j++) r[j] = b[j]-r[j];
460  beta1 = gmres_norm2(n,r);
461  eta = beta1/beta0;
462  };
463  };
464  iter++;
465  };
466 
467  if ( termcheck==CheckEachIter && MONITORLEVEL>1 )
468  fprintf(MONITOR," %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta);
469  if ( eta > tol && iter < max_iter ) goto restart;
470  else RCODE = ( eta <= tol ? 0 : 2 );
471  if ( MONITORLEVEL>1 && termcheck==CheckOnRestart )
472  fprintf(MONITOR," %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta);
473  data->mode = ( RCODE == 0 ? Solution : Final );
474  data->residuum = beta1; itlin_dataout(iter,n,y,data);
475 
476  errorexit:
477  if ( ERRORLEVEL > 0 && RCODE != 0 )
478  {
479  switch ( RCODE )
480  {
481  case 2:
482  fprintf(ERROR,"\n Error - Maximum allowed number of iterations exceeded\n");
483  break;
484  default :
485  fprintf(ERROR,"\n Error, code=%i\n",RCODE);
486  };
487  };
488  for (i=0;i<=i_max+1;i++) free(V[i]);
489  free(V); free(h); free(z); free(Hi); free(r); free(y0);
490 
491  info->precision = eta;
492  info->iter = iter;
493  info->nomatvec = nomatvec;
494  info->noprecl = nopreconl;
495  info->noprecr = nopreconr;
496 
497  // Forgotten:
498  free(q); free(w);
499  if (io_allocated) {
500  free(itlin_ioctl); itlin_ioctl=NULL;
501  };
502  if (data) free(data);
503 
504 }
505 
506 double gmres_norm2(int n, double *v)
507 { int i;
508  double rval = 0.0;
509  for (i=0;i<n;i++) rval += v[i]*v[i];
510  return sqrt( rval );
511 }
512 
513 int gmres_qrfact(int n, int lda, double *a, double *q, int ijob)
514 /*
515  Translation of Fortran routine DHEQR from file DGMRES of
516  library=slatec(slap) package.
517 */
518 { int i, j, j1, k,info=0;
519  double t, t1, t2, c, s;
520  if ( ijob <= 1 )
521  { for (k=0;k<n;k++)
522  { for (j=0;j<k-1;j++)
523  { i=2*j; c=q[i]; s=q[i+1];
524  j1=j+k*lda; t1=a[j1]; t2=a[j1+1];
525  a[j1]=c*t1-s*t2; a[j1+1]=s*t1+c*t2;
526  };
527  i=2*k; j1=k+lda*k; t1=a[j1]; t2=a[j1+1];
528  if ( t2==0.0 ) { c=1.0; s=0.0; }
529  else if ( fabs(t2)>=fabs(t1) )
530  { t=t1/t2; s=-1.0/sqrt(1.0+t*t); c=-s*t; }
531  else
532  { t=t2/t1; c=1.0/sqrt(1.0+t*t); s=-c*t; };
533  q[i]=c; q[i+1]=s; a[j1]=c*t1-s*t2;
534  if( a[j1]==0.0 ) info = k;
535  };
536  }
537  else
538  { for (k=0;k<n-1;k++)
539  {
540  i=2*k; j1=k+lda*(n-1);
541  t1=a[j1]; t2=a[j1+1];
542  c=q[i]; s=q[i+1];
543  a[j1]=c*t1-s*t2; a[j1+1]=s*t1 + c*t2;
544  };
545  info=0;
546  j1 = n-1+lda*(n-1); t1 = a[j1]; t2 = a[j1+1];
547  if ( t2==0.0 ) { c=1.0; s=0.0; }
548  else if ( fabs(t2)>=fabs(t1) )
549  { t=t1/t2; s=-1.0/sqrt(1.0+t*t); c=-s*t; }
550  else
551  { t=t2/t1; c=1.0/sqrt(1.0+t*t); s=-c*t; };
552  i=2*n-2; q[i]=c; q[i+1]=s; a[j1]=c*t1-s*t2;
553  if( a[j1]==0.0 ) info = n-1;
554  };
555  return info;
556 }
557 
558 void gmres_qrsolv(int n, int lda, double *a, double *q, double *b)
559 /*
560  Translation of Fortran routine DHELS from file DGMRES of
561  library=slatec(slap) package.
562 */
563 { int i, k, kb;
564  double t, t1, t2, c, s;
565  for (k=0;k<n;k++)
566  { i=2*k; c=q[i]; s=q[i+1];
567  t1=b[k]; t2=b[k+1];
568  b[k]=c*t1-s*t2; b[k+1]=s*t1+c*t2;
569  };
570  for (kb=0;kb<n;kb++)
571  {
572  k=n-1-kb; b[k] /= a[k+lda*k];
573  t = -b[k];
574  for (i=0;i<k;i++) b[i] += t*a[i+k*lda];
575  };
576  return;
577 }
static const double m
mass of electron in grams
void matvec(int n, double *x, double *b)
Used for gmres.