1 /* ./src_f77/ncgdrvd.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__9 = 9;
11 static integer c__1 = 1;
12 static integer c__3 = 3;
13 static integer c__0 = 0;
14 static integer c_n99 = -99;
15 static doublereal c_b15 = 0.;
16 static integer c__5 = 5;
17 static integer c_n1 = -1;
18 static doublereal c_b72 = 1.;
19 static doublereal c_b73 = -1.;
20
21 /* * /////////////////////////////////////////////////////////////////////////// */
22 /* * @file ncgdrvd.f */
23 /* * @author Michael Holst */
24 /* * @brief Driver for the nonlinear CG methods. */
25 /* * @version $Id: ncgdrvd.c,v 1.3 2010/08/12 05:45:24 fetk Exp $ */
26 /* * @attention */
27 /* * @verbatim */
28 /* * */
29 /* * PMG -- Parallel algebraic MultiGrid */
30 /* * Copyright (C) 1994-- Michael Holst. */
31 /* * */
32 /* * Michael Holst <mholst@math.ucsd.edu> */
33 /* * University of California, San Diego */
34 /* * Department of Mathematics, 5739 AP&M */
35 /* * 9500 Gilman Drive, Dept. 0112 */
36 /* * La Jolla, CA 92093-0112 USA */
37 /* * http://math.ucsd.edu/~mholst */
38 /* * */
39 /* * This file is part of PMG. */
40 /* * */
41 /* * PMG is free software; you can redistribute it and/or modify */
42 /* * it under the terms of the GNU General Public License as published by */
43 /* * the Free Software Foundation; either version 2 of the License, or */
44 /* * (at your option) any later version. */
45 /* * */
46 /* * PMG is distributed in the hope that it will be useful, */
47 /* * but WITHOUT ANY WARRANTY; without even the implied warranty of */
48 /* * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
49 /* * GNU General Public License for more details. */
50 /* * */
51 /* * You should have received a copy of the GNU General Public License */
52 /* * along with PMG; if not, write to the Free Software */
53 /* * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
54 /* * */
55 /* * Linking PMG statically or dynamically with other modules is making a */
56 /* * combined work based on PMG. Thus, the terms and conditions of the GNU */
57 /* * General Public License cover the whole combination. */
58 /* * */
59 /* * SPECIAL GPL EXCEPTION */
60 /* * In addition, as a special exception, the copyright holders of PMG */
61 /* * give you permission to combine the PMG program with free software */
62 /* * programs and libraries that are released under the GNU LGPL or with */
63 /* * code included in releases of ISIM, PMV, PyMOL, SMOL, VMD, and Vision. */
64 /* * Such combined software may be linked with PMG and redistributed together */
65 /* * in original or modified form as mere aggregation without requirement that */
66 /* * the entire work be under the scope of the GNU General Public License. */
67 /* * This special exception permission is also extended to any software listed */
68 /* * in the SPECIAL GPL EXCEPTION clauses by the FEtk and APBS libraries. */
69 /* * */
70 /* * Note that people who make modified versions of PMG are not obligated */
71 /* * to grant this special exception for their modified versions; it is */
72 /* * their choice whether to do so. The GNU General Public License gives */
73 /* * permission to release a modified version without this exception; this */
74 /* * exception also makes it possible to release a modified version which */
75 /* * carries forward this exception. */
76 /* * */
77 /* * @endverbatim */
78 /* * /////////////////////////////////////////////////////////////////////////// */
ncghsdriv_(integer * iparm,doublereal * rparm,integer * iwork,doublereal * rwork,doublereal * u,doublereal * xf,doublereal * yf,doublereal * zf,doublereal * gxcf,doublereal * gycf,doublereal * gzcf,doublereal * a1cf,doublereal * a2cf,doublereal * a3cf,doublereal * ccf,doublereal * fcf,doublereal * tcf)79 /* Subroutine */ int ncghsdriv_(integer *iparm, doublereal *rparm, integer *
80 iwork, doublereal *rwork, doublereal *u, doublereal *xf, doublereal *
81 yf, doublereal *zf, doublereal *gxcf, doublereal *gycf, doublereal *
82 gzcf, doublereal *a1cf, doublereal *a2cf, doublereal *a3cf,
83 doublereal *ccf, doublereal *fcf, doublereal *tcf)
84 {
85 /* Builtin functions */
86 integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
87 e_wsle(void);
88
89 /* Local variables */
90 static integer n;
91 extern /* Subroutine */ int ncghsdriv2_(integer *, doublereal *, integer *
92 , integer *, integer *, doublereal *, integer *, doublereal *,
93 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
94 doublereal *, doublereal *, doublereal *, doublereal *,
95 doublereal *, doublereal *, doublereal *, doublereal *,
96 doublereal *, doublereal *, doublereal *, doublereal *);
97 static integer nx, ny, nz, k_w0__, k_ac__, k_cc__, k_fc__, k_iz__, n_iz__,
98 nlev, niwk, nrwk, k_ipc__, n_ipc__, k_rpc__, n_rpc__, ierror,
99 iretot, iintot;
100
101 /* Fortran I/O blocks */
102 static cilist io___13 = { 0, 6, 0, 0, 0 };
103 static cilist io___14 = { 0, 6, 0, 0, 0 };
104
105
106 /* * ********************************************************************* */
107 /* * purpose: */
108 /* * */
109 /* * linear/nonlinear conjugate gradient driver (fletcher-reeves). */
110 /* * */
111 /* * author: michael holst */
112 /* * ********************************************************************* */
113 /* * */
114 /* * *** input parameters *** */
115 /* * */
116 /* * *** variables returned from mgsz *** */
117 /* * */
118 /* * *** misc variables *** */
119 /* * */
120 /* * *** decode some parameters *** */
121 /* Parameter adjustments */
122 --tcf;
123 --fcf;
124 --ccf;
125 --a3cf;
126 --a2cf;
127 --a1cf;
128 --gzcf;
129 --gycf;
130 --gxcf;
131 --zf;
132 --yf;
133 --xf;
134 --u;
135 --rwork;
136 --iwork;
137 --rparm;
138 --iparm;
139
140 /* Function Body */
141 nrwk = iparm[1];
142 niwk = iparm[2];
143 nx = iparm[3];
144 ny = iparm[4];
145 nz = iparm[5];
146 nlev = iparm[6];
147 n = nx * ny * nz;
148 n_iz__ = (nlev + 1) * 10;
149 n_ipc__ = (nlev + 1) * 100;
150 n_rpc__ = (nlev + 1) * 100;
151 /* * */
152 /* * *** compute required work array sizes *** */
153 iintot = n_iz__ + n_ipc__;
154 iretot = n_rpc__ + n * 3 + (n << 2);
155 /* * */
156 /* * *** some more checks on input *** */
157 if (nrwk < iretot || niwk < iintot) {
158 s_wsle(&io___13);
159 do_lio(&c__9, &c__1, "% NCGHSDRIV: real work space must be: ", (
160 ftnlen)41);
161 do_lio(&c__3, &c__1, (char *)&iretot, (ftnlen)sizeof(integer));
162 e_wsle();
163 s_wsle(&io___14);
164 do_lio(&c__9, &c__1, "% NCGHSDRIV: integer work space must be: ", (
165 ftnlen)41);
166 do_lio(&c__3, &c__1, (char *)&iintot, (ftnlen)sizeof(integer));
167 e_wsle();
168 ierror = -3;
169 iparm[51] = ierror;
170 return 0;
171 }
172 /* * */
173 /* * *** iwork offsets *** */
174 k_iz__ = 1;
175 k_ipc__ = k_iz__ + n_iz__;
176 /* * */
177 /* * *** rwork offsets *** */
178 k_rpc__ = 1;
179 k_cc__ = k_rpc__ + n_rpc__;
180 k_fc__ = k_cc__ + n;
181 k_w0__ = k_fc__ + n;
182 k_ac__ = k_w0__ + n;
183 /* * ***k_ac_after = k_ac + 4*n */
184 /* * */
185 /* * *** call the multigrid driver *** */
186 ncghsdriv2_(&iparm[1], &rparm[1], &nx, &ny, &nz, &u[1], &iwork[k_iz__], &
187 rwork[k_w0__], &iwork[k_ipc__], &rwork[k_rpc__], &rwork[k_ac__], &
188 rwork[k_cc__], &rwork[k_fc__], &xf[1], &yf[1], &zf[1], &gxcf[1], &
189 gycf[1], &gzcf[1], &a1cf[1], &a2cf[1], &a3cf[1], &ccf[1], &fcf[1],
190 &tcf[1]);
191 /* * */
192 /* * *** return and end *** */
193 return 0;
194 } /* ncghsdriv_ */
195
ncghsdriv2_(integer * iparm,doublereal * rparm,integer * nx,integer * ny,integer * nz,doublereal * u,integer * iz,doublereal * w0,integer * ipc,doublereal * rpc,doublereal * ac,doublereal * cc,doublereal * fc,doublereal * xf,doublereal * yf,doublereal * zf,doublereal * gxcf,doublereal * gycf,doublereal * gzcf,doublereal * a1cf,doublereal * a2cf,doublereal * a3cf,doublereal * ccf,doublereal * fcf,doublereal * tcf)196 /* Subroutine */ int ncghsdriv2_(integer *iparm, doublereal *rparm, integer *
197 nx, integer *ny, integer *nz, doublereal *u, integer *iz, doublereal *
198 w0, integer *ipc, doublereal *rpc, doublereal *ac, doublereal *cc,
199 doublereal *fc, doublereal *xf, doublereal *yf, doublereal *zf,
200 doublereal *gxcf, doublereal *gycf, doublereal *gzcf, doublereal *
201 a1cf, doublereal *a2cf, doublereal *a3cf, doublereal *ccf, doublereal
202 *fcf, doublereal *tcf)
203 {
204 /* Format strings */
205 static char fmt_100[] = "(\002 ];\002,4(\002 \002,a7,\002=\002,1pe9.3"
206 ",\002;\002))";
207
208 /* System generated locals */
209 doublereal d__1;
210
211 /* Builtin functions */
212 integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
213 e_wsle(void), s_wsfe(cilist *), do_fio(integer *, char *, ftnlen),
214 e_wsfe(void);
215
216 /* Local variables */
217 extern /* Subroutine */ int buildops_(integer *, integer *, integer *,
218 integer *, integer *, integer *, integer *, integer *, integer *,
219 integer *, integer *, integer *, integer *, doublereal *,
220 doublereal *, doublereal *, doublereal *, doublereal *,
221 doublereal *, doublereal *, doublereal *, doublereal *,
222 doublereal *, doublereal *, doublereal *, doublereal *,
223 doublereal *, doublereal *, doublereal *, doublereal *),
224 buildstr_(integer *, integer *, integer *, integer *, integer *);
225 static doublereal bf, oh;
226 static integer nu1, nu2, ido, iok, mode, ilev, nlev, iinfo, mgkey, itmax,
227 ipkey, iters, istop;
228 extern /* Subroutine */ int tstop_(doublereal *, doublereal *, doublereal
229 *);
230 static doublereal omegal;
231 static integer mgdisc;
232 static doublereal omegan;
233 static integer mgcoar;
234 extern doublereal epsmac_(integer *);
235 extern /* Subroutine */ int cghsgo_(integer *, integer *, integer *,
236 doublereal *, doublereal *, doublereal *, doublereal *,
237 doublereal *, doublereal *, doublereal *, integer *, integer *,
238 integer *, integer *, integer *, integer *, doublereal *,
239 doublereal *, doublereal *, integer *, doublereal *, doublereal *,
240 doublereal *, doublereal *, doublereal *), fbound_(integer *,
241 integer *, integer *, integer *, doublereal *, doublereal *,
242 doublereal *, doublereal *), matvec_(integer *, integer *,
243 integer *, integer *, doublereal *, doublereal *, doublereal *,
244 doublereal *, doublereal *);
245 static integer ibound;
246 static doublereal epsiln;
247 static integer mgprol, ierror, mgsolv;
248 static doublereal errtol, tsolve;
249 extern /* Subroutine */ int tstart_(doublereal *, doublereal *), prtstp_(
250 integer *, integer *, doublereal *, doublereal *, doublereal *),
251 fbound00_(integer *, integer *, integer *, doublereal *);
252 static doublereal pc_dumm__;
253 extern /* Subroutine */ int ncghsgo_(integer *, integer *, integer *,
254 doublereal *, doublereal *, doublereal *, doublereal *,
255 doublereal *, doublereal *, doublereal *, integer *, integer *,
256 integer *, integer *, integer *, integer *, doublereal *,
257 doublereal *, doublereal *, integer *, doublereal *, doublereal *,
258 doublereal *, doublereal *, doublereal *), nmatvec_(integer *,
259 integer *, integer *, integer *, doublereal *, doublereal *,
260 doublereal *, doublereal *, doublereal *, doublereal *);
261 static doublereal tsetupc, tsetupf;
262
263 /* Fortran I/O blocks */
264 static cilist io___44 = { 0, 6, 0, 0, 0 };
265 static cilist io___47 = { 0, 6, 0, 0, 0 };
266 static cilist io___48 = { 0, 6, 0, 0, 0 };
267 static cilist io___53 = { 0, 6, 0, 0, 0 };
268 static cilist io___55 = { 0, 6, 0, 0, 0 };
269 static cilist io___56 = { 0, 6, 0, fmt_100, 0 };
270
271
272 /* * ********************************************************************* */
273 /* * purpose: */
274 /* * */
275 /* * linear/nonlinear conjugate gradient driver (fletcher-reeves). */
276 /* * */
277 /* * author: michael holst */
278 /* * ********************************************************************* */
279 /* * */
280 /* * *** input parameters *** */
281 /* * */
282 /* * *** misc variables *** */
283 /* * */
284 /* * *** decode the iparm array *** */
285 /* Parameter adjustments */
286 --tcf;
287 --fcf;
288 --ccf;
289 --a3cf;
290 --a2cf;
291 --a1cf;
292 --gzcf;
293 --gycf;
294 --gxcf;
295 --zf;
296 --yf;
297 --xf;
298 --fc;
299 --cc;
300 --ac;
301 --rpc;
302 --ipc;
303 --w0;
304 --iz;
305 --u;
306 --rparm;
307 --iparm;
308
309 /* Function Body */
310 nlev = iparm[6];
311 nu1 = iparm[7];
312 nu2 = iparm[8];
313 mgkey = iparm[9];
314 itmax = iparm[10];
315 istop = iparm[11];
316 iinfo = iparm[12];
317 ipkey = iparm[14];
318 mode = iparm[16];
319 mgprol = iparm[17];
320 mgcoar = iparm[18];
321 mgdisc = iparm[19];
322 mgsolv = iparm[21];
323 errtol = rparm[1];
324 omegal = rparm[9];
325 omegan = rparm[10];
326 /* * */
327 /* * *** intitialize the iteration timer *** */
328 prtstp_(&c__0, &c_n99, &c_b15, &c_b15, &c_b15);
329 /* * */
330 /* * *** build the multigrid data structure in iz *** */
331 buildstr_(nx, ny, nz, &nlev, &iz[1]);
332 /* * */
333 /* * *** start timer *** */
334 tstart_(&bf, &oh);
335 /* * */
336 /* * *** build op and rhs on fine grid *** */
337 ido = 0;
338 buildops_(nx, ny, nz, &nlev, &ipkey, &iinfo, &ido, &iz[1], &mgprol, &
339 mgcoar, &mgsolv, &mgdisc, &ipc[1], &rpc[1], &pc_dumm__, &ac[1], &
340 cc[1], &fc[1], &xf[1], &yf[1], &zf[1], &gxcf[1], &gycf[1], &gzcf[
341 1], &a1cf[1], &a2cf[1], &a3cf[1], &ccf[1], &fcf[1], &tcf[1]);
342 /* * */
343 /* * *** stop timer *** */
344 tstop_(&bf, &oh, &tsetupf);
345 s_wsle(&io___44);
346 do_lio(&c__9, &c__1, "% NCGHSDRIV2: fine problem setup time: ", (ftnlen)
347 39);
348 do_lio(&c__5, &c__1, (char *)&tsetupf, (ftnlen)sizeof(doublereal));
349 e_wsle();
350 tsetupc = 0.;
351 /* * */
352 /* * ****************************************************************** */
353 /* * *** this overwrites the rhs array provided by pde specification */
354 /* * ****** compute an algebraically produced rhs for the given tcf *** */
355 if (istop == 4 || istop == 5) {
356 if (mode == 1 || mode == 2) {
357 nmatvec_(nx, ny, nz, &ipc[1], &rpc[1], &ac[1], &cc[1], &tcf[1], &
358 fc[1], &w0[1]);
359 } else {
360 matvec_(nx, ny, nz, &ipc[1], &rpc[1], &ac[1], &cc[1], &tcf[1], &
361 fc[1]);
362 }
363 }
364 /* * ****************************************************************** */
365 /* * */
366 /* * *** determine machine epsilon *** */
367 epsiln = epsmac_(&c__0);
368 /* * */
369 /* * *** impose zero dirichlet boundary conditions (now in source fcn) *** */
370 fbound00_(nx, ny, nz, &u[1]);
371 /* * */
372 /* * *** MATLAB *** */
373 s_wsle(&io___47);
374 do_lio(&c__9, &c__1, " cg = [ ", (ftnlen)8);
375 e_wsle();
376 /* * */
377 /* * *** start timer *** */
378 tstart_(&bf, &oh);
379 /* * */
380 /* * *** call specified multigrid method *** */
381 if (mode == 0 || mode == 2) {
382 s_wsle(&io___48);
383 do_lio(&c__9, &c__1, "% NCGHSDRIV2: linear mode...", (ftnlen)28);
384 e_wsle();
385 iok = 1;
386 ilev = 1;
387 cghsgo_(nx, ny, nz, &u[1], &w0[1], &a1cf[1], &a2cf[1], &a3cf[1], &ccf[
388 1], &fcf[1], &istop, &itmax, &iters, &ierror, &iok, &iinfo, &
389 epsiln, &errtol, &omegal, &ipc[1], &rpc[1], &ac[1], &cc[1], &
390 fc[1], &tcf[1]);
391 }
392 if (mode == 1 || mode == 2) {
393 s_wsle(&io___53);
394 do_lio(&c__9, &c__1, "% NCGHSDRIV2: nonlinear mode...", (ftnlen)31);
395 e_wsle();
396 iok = 1;
397 ilev = 1;
398 ncghsgo_(nx, ny, nz, &u[1], &w0[1], &a1cf[1], &a2cf[1], &a3cf[1], &
399 ccf[1], &fcf[1], &istop, &itmax, &iters, &ierror, &iok, &
400 iinfo, &epsiln, &errtol, &omegan, &ipc[1], &rpc[1], &ac[1], &
401 cc[1], &fc[1], &tcf[1]);
402 }
403 /* * */
404 /* * *** stop timer *** */
405 tstop_(&bf, &oh, &tsolve);
406 s_wsle(&io___55);
407 do_lio(&c__9, &c__1, "% NCGHSDRIV2: solve time: ", (ftnlen)26);
408 do_lio(&c__5, &c__1, (char *)&tsolve, (ftnlen)sizeof(doublereal));
409 e_wsle();
410 /* * */
411 /* * *** MATLAB *** */
412 s_wsfe(&io___56);
413 do_fio(&c__1, "cg_sf", (ftnlen)5);
414 do_fio(&c__1, (char *)&tsetupf, (ftnlen)sizeof(doublereal));
415 do_fio(&c__1, "cg_sc", (ftnlen)5);
416 do_fio(&c__1, (char *)&tsetupc, (ftnlen)sizeof(doublereal));
417 do_fio(&c__1, "cg_st", (ftnlen)5);
418 d__1 = tsetupf + tsetupc;
419 do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal));
420 do_fio(&c__1, "cg_so", (ftnlen)5);
421 do_fio(&c__1, (char *)&tsolve, (ftnlen)sizeof(doublereal));
422 e_wsfe();
423 /* * */
424 /* * *** restore boundary conditions *** */
425 ibound = 1;
426 fbound_(&ibound, nx, ny, nz, &u[1], &gxcf[1], &gycf[1], &gzcf[1]);
427 /* * */
428 /* * *** return and end *** */
429 return 0;
430 } /* ncghsdriv2_ */
431
ncghsgo_(integer * nx,integer * ny,integer * nz,doublereal * x,doublereal * r__,doublereal * p,doublereal * ap,doublereal * zk,doublereal * zkp1,doublereal * tmp,integer * istop,integer * itmax,integer * iters,integer * ierror,integer * iok,integer * iinfo,doublereal * epsiln,doublereal * errtol,doublereal * omega,integer * ipc,doublereal * rpc,doublereal * ac,doublereal * cc,doublereal * fc,doublereal * tru)432 /* Subroutine */ int ncghsgo_(integer *nx, integer *ny, integer *nz,
433 doublereal *x, doublereal *r__, doublereal *p, doublereal *ap,
434 doublereal *zk, doublereal *zkp1, doublereal *tmp, integer *istop,
435 integer *itmax, integer *iters, integer *ierror, integer *iok,
436 integer *iinfo, doublereal *epsiln, doublereal *errtol, doublereal *
437 omega, integer *ipc, doublereal *rpc, doublereal *ac, doublereal *cc,
438 doublereal *fc, doublereal *tru)
439 {
440 /* Format strings */
441 static char fmt_100[] = "(a,(2x,\002 [\002,i3,\002,\002,i3,\002,\002,i3"
442 ",\002] \002))";
443
444 /* System generated locals */
445 doublereal d__1;
446
447 /* Builtin functions */
448 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
449 double sqrt(doublereal);
450 integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
451 e_wsle(void);
452
453 /* Local variables */
454 extern /* Subroutine */ int linesearch_(integer *, integer *, integer *,
455 doublereal *, integer *, doublereal *, doublereal *, doublereal *,
456 doublereal *, doublereal *, doublereal *, doublereal *,
457 doublereal *, doublereal *, doublereal *);
458 static doublereal beta;
459 extern doublereal xdot_(integer *, integer *, integer *, doublereal *,
460 doublereal *);
461 static doublereal rhok1, rhok2;
462 extern doublereal xnrm1_(integer *, integer *, integer *, doublereal *),
463 xnrm2_(integer *, integer *, integer *, doublereal *);
464 static doublereal alpha;
465 extern /* Subroutine */ int xscal_(integer *, integer *, integer *,
466 doublereal *, doublereal *);
467 static doublereal rsden, rsnrm;
468 extern /* Subroutine */ int xcopy_(integer *, integer *, integer *,
469 doublereal *, doublereal *), xaxpy_(integer *, integer *, integer
470 *, doublereal *, doublereal *, doublereal *), azeros_(integer *,
471 integer *, integer *, doublereal *), prtini_(integer *);
472 static doublereal orsnrm;
473 extern /* Subroutine */ int prtstp_(integer *, integer *, doublereal *,
474 doublereal *, doublereal *), nmatvec_(integer *, integer *,
475 integer *, integer *, doublereal *, doublereal *, doublereal *,
476 doublereal *, doublereal *, doublereal *), nmresid_(integer *,
477 integer *, integer *, integer *, doublereal *, doublereal *,
478 doublereal *, doublereal *, doublereal *, doublereal *,
479 doublereal *);
480
481 /* Fortran I/O blocks */
482 static cilist io___58 = { 0, 6, 0, fmt_100, 0 };
483 static cilist io___60 = { 0, 6, 0, 0, 0 };
484 static cilist io___61 = { 0, 6, 0, 0, 0 };
485 static cilist io___68 = { 0, 6, 0, 0, 0 };
486
487
488 /* * ********************************************************************* */
489 /* * purpose: */
490 /* * */
491 /* * nonlinear conjugate gradients (fletcher-reeves). */
492 /* * */
493 /* * author: michael holst */
494 /* * ********************************************************************* */
495 /* * */
496 /* * *** other declarations *** */
497 /* * */
498 /* mdir 0 0 */
499 /* * */
500 /* * *** do some i/o if requested *** */
501 /* Parameter adjustments */
502 --tru;
503 --fc;
504 --cc;
505 --ac;
506 --rpc;
507 --ipc;
508 --tmp;
509 --zkp1;
510 --zk;
511 --ap;
512 --p;
513 --r__;
514 --x;
515
516 /* Function Body */
517 if (*iinfo != 0) {
518 s_wsfe(&io___58);
519 do_fio(&c__1, "% NCGHSGO: starting:", (ftnlen)20);
520 do_fio(&c__1, (char *)&(*nx), (ftnlen)sizeof(integer));
521 do_fio(&c__1, (char *)&(*ny), (ftnlen)sizeof(integer));
522 do_fio(&c__1, (char *)&(*nz), (ftnlen)sizeof(integer));
523 e_wsfe();
524 }
525 /* * */
526 /* * *** initial wall clock *** */
527 prtini_(istop);
528 prtstp_(iok, &c_n1, &c_b15, &c_b15, &c_b15);
529 /* * */
530 /* * ************************************************************** */
531 /* * *** note: if (iok.ne.0) then: use a stopping test. *** */
532 /* * *** else: use just the itmax to stop iteration. *** */
533 /* * ************************************************************** */
534 /* * *** istop=0 most efficient (whatever it is) *** */
535 /* * *** istop=1 relative residual *** */
536 /* * *** istop=2 rms difference of successive iterates *** */
537 /* * *** istop=3 relative true error (provided for testing) *** */
538 /* * ************************************************************** */
539 /* * */
540 /* * *** compute denominator for stopping criterion *** */
541 if (*istop == 0) {
542 rsden = 1.;
543 } else if (*istop == 1) {
544 /* * *** compute initial residual with zero initial guess *** */
545 /* * *** this is analogous to the linear case where one can *** */
546 /* * *** simply take norm of rhs for a zero initial guess *** */
547 azeros_(nx, ny, nz, &tmp[1]);
548 nmresid_(nx, ny, nz, &ipc[1], &rpc[1], &ac[1], &cc[1], &fc[1], &tmp[1]
549 , &r__[1], &zk[1]);
550 rsden = xnrm1_(nx, ny, nz, &r__[1]);
551 } else if (*istop == 2) {
552 rsden = sqrt((doublereal) ((*nx - 2) * (*ny - 2) * (*nz - 2)));
553 } else if (*istop == 3) {
554 rsden = xnrm2_(nx, ny, nz, &tru[1]);
555 } else if (*istop == 4) {
556 rsden = xnrm2_(nx, ny, nz, &tru[1]);
557 } else if (*istop == 5) {
558 nmatvec_(nx, ny, nz, &ipc[1], &rpc[1], &ac[1], &cc[1], &tru[1], &r__[
559 1], &zk[1]);
560 rsden = sqrt(xdot_(nx, ny, nz, &tru[1], &r__[1]));
561 } else {
562 s_wsle(&io___60);
563 do_lio(&c__9, &c__1, "% NCGHSGO: bad istop value... ", (ftnlen)30);
564 e_wsle();
565 }
566 if (rsden == 0.) {
567 rsden = 1.;
568 s_wsle(&io___61);
569 do_lio(&c__9, &c__1, "% NCGHSGO: rhs is zero ", (ftnlen)23);
570 e_wsle();
571 }
572 rsnrm = rsden;
573 orsnrm = rsnrm;
574 prtstp_(iok, &c__0, &rsnrm, &rsden, &orsnrm);
575 /* * */
576 /* * *** now compute residual with the initial guess *** */
577 nmresid_(nx, ny, nz, &ipc[1], &rpc[1], &ac[1], &cc[1], &fc[1], &x[1], &
578 r__[1], &zk[1]);
579 /* * */
580 /* * *** setup for the looping *** */
581 *iters = 0;
582 L30:
583 /* * */
584 /* * *** save iterate if stop test will use it on next iter *** */
585 if (*istop == 2) {
586 xcopy_(nx, ny, nz, &x[1], &tru[1]);
587 }
588 /* * */
589 /* * *** form new direction vector from old one and residual *** */
590 rhok2 = xdot_(nx, ny, nz, &r__[1], &r__[1]);
591 if (*iters == 0) {
592 xcopy_(nx, ny, nz, &r__[1], &p[1]);
593 } else {
594 beta = rhok2 / rhok1;
595 d__1 = 1. / beta;
596 xaxpy_(nx, ny, nz, &d__1, &r__[1], &p[1]);
597 xscal_(nx, ny, nz, &beta, &p[1]);
598 }
599 /* * */
600 /* * *** nonlinear case: do a line search *** */
601 /* * *** (note: "ap,zk,zkp1" passed back from line search as desired) *** */
602 xcopy_(nx, ny, nz, &r__[1], &tmp[1]);
603 linesearch_(nx, ny, nz, &alpha, &ipc[1], &rpc[1], &ac[1], &cc[1], &fc[1],
604 &p[1], &x[1], &tmp[1], &ap[1], &zk[1], &zkp1[1]);
605 /* * */
606 /* * *** save rhok2 for next iteration *** */
607 rhok1 = rhok2;
608 /* * */
609 /* * *** update solution in direction p of length alpha *** */
610 xaxpy_(nx, ny, nz, &alpha, &p[1], &x[1]);
611 /* * */
612 /* * *** update residual *** */
613 d__1 = -alpha;
614 xaxpy_(nx, ny, nz, &d__1, &ap[1], &r__[1]);
615 xaxpy_(nx, ny, nz, &c_b72, &zk[1], &r__[1]);
616 xaxpy_(nx, ny, nz, &c_b73, &zkp1[1], &r__[1]);
617 /* * */
618 /* * ***** *** switch to descent if necessary *** */
619 /* * ***** itmax_d = 10 */
620 /* * ***** iter_d = 0 */
621 /* * ***** alpha = -1.0d0 */
622 /* * ***** rsnrm_tmp = xnrm1(nx,ny,nz,r) */
623 /* * ***** 18 continue */
624 /* * ***** if ((rsnrm_tmp.lt.rsnrm).or.(iter_d.gt.itmax_d)) then */
625 /* * ***** print*,'% finished with descent: r_o,r_n',rsnrm,rsnrm_tmp */
626 /* * ***** if (iter_d .gt. 0) call xcopy(nx,ny,nz,tmp4,r) */
627 /* * ***** goto 19 */
628 /* * ***** endif */
629 /* * ***** print*,'% trying a descent: r_o,r_n',rsnrm,rsnrm_tmp */
630 /* * ***** call xcopy(nx,ny,nz,tmp2,tmp) */
631 /* * ***** call xaxpy(nx,ny,nz,alpha,tmp3,tmp) */
632 /* * ***** call nmresid(nx,ny,nz,ipc,rpc,ac,cc,fc,tmp,tmp4,ap) */
633 /* * ***** rsnrm_tmp = xnrm1(nx,ny,nz,tmp4) */
634 /* * ***** alpha = alpha / 2.0d0 */
635 /* * ***** iter_d = iter_d + 1 */
636 /* * ***** goto 18 */
637 /* * ***** 19 continue */
638 /* * */
639 /* * *** some bookkeeping *** */
640 ++(*iters);
641 /* * */
642 /* * *** compute/check the current stopping test *** */
643 orsnrm = rsnrm;
644 if (*istop == 0) {
645 rsnrm = xnrm1_(nx, ny, nz, &r__[1]);
646 } else if (*istop == 1) {
647 rsnrm = xnrm1_(nx, ny, nz, &r__[1]);
648 } else if (*istop == 2) {
649 xcopy_(nx, ny, nz, &tru[1], &tmp[1]);
650 xaxpy_(nx, ny, nz, &c_b73, &x[1], &tmp[1]);
651 rsnrm = xnrm1_(nx, ny, nz, &tmp[1]);
652 } else if (*istop == 3) {
653 xcopy_(nx, ny, nz, &tru[1], &tmp[1]);
654 xaxpy_(nx, ny, nz, &c_b73, &x[1], &tmp[1]);
655 rsnrm = xnrm2_(nx, ny, nz, &tmp[1]);
656 } else if (*istop == 4) {
657 xcopy_(nx, ny, nz, &tru[1], &tmp[1]);
658 xaxpy_(nx, ny, nz, &c_b73, &x[1], &tmp[1]);
659 rsnrm = xnrm2_(nx, ny, nz, &tmp[1]);
660 } else if (*istop == 5) {
661 xcopy_(nx, ny, nz, &tru[1], &tmp[1]);
662 xaxpy_(nx, ny, nz, &c_b73, &x[1], &tmp[1]);
663 nmatvec_(nx, ny, nz, &ipc[1], &rpc[1], &ac[1], &cc[1], &tmp[1], &zk[1]
664 , &zkp1[1]);
665 rsnrm = sqrt(xdot_(nx, ny, nz, &tmp[1], &zk[1]));
666 } else {
667 s_wsle(&io___68);
668 do_lio(&c__9, &c__1, "% NCGHSGO: bad istop value... ", (ftnlen)30);
669 e_wsle();
670 }
671 prtstp_(iok, iters, &rsnrm, &rsden, &orsnrm);
672 if (rsnrm / rsden <= *errtol) {
673 goto L99;
674 }
675 if (*iters >= *itmax) {
676 goto L99;
677 }
678 goto L30;
679 /* * */
680 /* * *** return and end *** */
681 L99:
682 return 0;
683 } /* ncghsgo_ */
684
cghsgo_(integer * nx,integer * ny,integer * nz,doublereal * x,doublereal * r__,doublereal * p,doublereal * ap,doublereal * zk,doublereal * zkp1,doublereal * tmp,integer * istop,integer * itmax,integer * iters,integer * ierror,integer * iok,integer * iinfo,doublereal * epsiln,doublereal * errtol,doublereal * omega,integer * ipc,doublereal * rpc,doublereal * ac,doublereal * cc,doublereal * fc,doublereal * tru)685 /* Subroutine */ int cghsgo_(integer *nx, integer *ny, integer *nz,
686 doublereal *x, doublereal *r__, doublereal *p, doublereal *ap,
687 doublereal *zk, doublereal *zkp1, doublereal *tmp, integer *istop,
688 integer *itmax, integer *iters, integer *ierror, integer *iok,
689 integer *iinfo, doublereal *epsiln, doublereal *errtol, doublereal *
690 omega, integer *ipc, doublereal *rpc, doublereal *ac, doublereal *cc,
691 doublereal *fc, doublereal *tru)
692 {
693 /* Format strings */
694 static char fmt_100[] = "(a,(2x,\002 [\002,i3,\002,\002,i3,\002,\002,i3"
695 ",\002] \002))";
696
697 /* System generated locals */
698 doublereal d__1;
699
700 /* Builtin functions */
701 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
702 double sqrt(doublereal);
703 integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
704 e_wsle(void);
705
706 /* Local variables */
707 static doublereal pap, beta;
708 extern doublereal xdot_(integer *, integer *, integer *, doublereal *,
709 doublereal *);
710 static doublereal rhok1, rhok2;
711 extern doublereal xnrm1_(integer *, integer *, integer *, doublereal *),
712 xnrm2_(integer *, integer *, integer *, doublereal *);
713 static doublereal alpha;
714 extern /* Subroutine */ int xscal_(integer *, integer *, integer *,
715 doublereal *, doublereal *);
716 static doublereal rsden, rsnrm;
717 extern /* Subroutine */ int xcopy_(integer *, integer *, integer *,
718 doublereal *, doublereal *), xaxpy_(integer *, integer *, integer
719 *, doublereal *, doublereal *, doublereal *), matvec_(integer *,
720 integer *, integer *, integer *, doublereal *, doublereal *,
721 doublereal *, doublereal *, doublereal *), mresid_(integer *,
722 integer *, integer *, integer *, doublereal *, doublereal *,
723 doublereal *, doublereal *, doublereal *, doublereal *), prtini_(
724 integer *);
725 static doublereal orsnrm;
726 extern /* Subroutine */ int prtstp_(integer *, integer *, doublereal *,
727 doublereal *, doublereal *);
728
729 /* Fortran I/O blocks */
730 static cilist io___69 = { 0, 6, 0, fmt_100, 0 };
731 static cilist io___71 = { 0, 6, 0, 0, 0 };
732 static cilist io___72 = { 0, 6, 0, 0, 0 };
733 static cilist io___80 = { 0, 6, 0, 0, 0 };
734
735
736 /* * ********************************************************************* */
737 /* * purpose: */
738 /* * */
739 /* * linear conjugate gradients (hestenes-steifel). */
740 /* * */
741 /* * author: michael holst */
742 /* * ********************************************************************* */
743 /* * */
744 /* * *** other declarations *** */
745 /* * */
746 /* mdir 0 0 */
747 /* * */
748 /* * *** do some i/o if requested *** */
749 /* Parameter adjustments */
750 --tru;
751 --fc;
752 --cc;
753 --ac;
754 --rpc;
755 --ipc;
756 --tmp;
757 --zkp1;
758 --zk;
759 --ap;
760 --p;
761 --r__;
762 --x;
763
764 /* Function Body */
765 if (*iinfo != 0) {
766 s_wsfe(&io___69);
767 do_fio(&c__1, "% CGHSGO: starting: ", (ftnlen)20);
768 do_fio(&c__1, (char *)&(*nx), (ftnlen)sizeof(integer));
769 do_fio(&c__1, (char *)&(*ny), (ftnlen)sizeof(integer));
770 do_fio(&c__1, (char *)&(*nz), (ftnlen)sizeof(integer));
771 e_wsfe();
772 }
773 /* * */
774 /* * *** initial wall clock *** */
775 prtini_(istop);
776 prtstp_(iok, &c_n1, &c_b15, &c_b15, &c_b15);
777 /* * */
778 /* * ************************************************************** */
779 /* * *** note: if (iok.ne.0) then: use a stopping test. *** */
780 /* * *** else: use just the itmax to stop iteration. *** */
781 /* * ************************************************************** */
782 /* * *** istop=0 most efficient (whatever it is) *** */
783 /* * *** istop=1 relative residual *** */
784 /* * *** istop=2 rms difference of successive iterates *** */
785 /* * *** istop=3 relative true error (provided for testing) *** */
786 /* * ************************************************************** */
787 /* * */
788 /* * *** compute denominator for stopping criterion *** */
789 if (*istop == 0) {
790 rsden = 1.;
791 } else if (*istop == 1) {
792 rsden = xnrm1_(nx, ny, nz, &fc[1]);
793 } else if (*istop == 2) {
794 rsden = sqrt((doublereal) ((*nx - 2) * (*ny - 2) * (*nz - 2)));
795 } else if (*istop == 3) {
796 rsden = xnrm2_(nx, ny, nz, &tru[1]);
797 } else if (*istop == 4) {
798 rsden = xnrm2_(nx, ny, nz, &tru[1]);
799 } else if (*istop == 5) {
800 matvec_(nx, ny, nz, &ipc[1], &rpc[1], &ac[1], &cc[1], &tru[1], &r__[1]
801 );
802 rsden = sqrt(xdot_(nx, ny, nz, &tru[1], &r__[1]));
803 } else {
804 s_wsle(&io___71);
805 do_lio(&c__9, &c__1, "% CGHSGO: bad istop value... ", (ftnlen)29);
806 e_wsle();
807 }
808 if (rsden == 0.) {
809 rsden = 1.;
810 s_wsle(&io___72);
811 do_lio(&c__9, &c__1, "% CGHSGO: rhs is zero ", (ftnlen)22);
812 e_wsle();
813 }
814 rsnrm = rsden;
815 orsnrm = rsnrm;
816 prtstp_(iok, &c__0, &rsnrm, &rsden, &orsnrm);
817 /* * */
818 /* * *** now compute residual with the initial guess *** */
819 mresid_(nx, ny, nz, &ipc[1], &rpc[1], &ac[1], &cc[1], &fc[1], &x[1], &r__[
820 1]);
821 /* * */
822 /* * *** setup for the looping *** */
823 *iters = 0;
824 L30:
825 /* * */
826 /* * *** save iterate if stop test will use it on next iter *** */
827 if (*istop == 2) {
828 xcopy_(nx, ny, nz, &x[1], &tru[1]);
829 }
830 /* * */
831 /* * *** form new direction vector from old one and residual *** */
832 rhok2 = xdot_(nx, ny, nz, &r__[1], &r__[1]);
833 if (*iters == 0) {
834 xcopy_(nx, ny, nz, &r__[1], &p[1]);
835 } else {
836 beta = rhok2 / rhok1;
837 d__1 = 1. / beta;
838 xaxpy_(nx, ny, nz, &d__1, &r__[1], &p[1]);
839 xscal_(nx, ny, nz, &beta, &p[1]);
840 }
841 /* * */
842 /* * *** linear case: alpha which minimizes energy norm of error *** */
843 matvec_(nx, ny, nz, &ipc[1], &rpc[1], &ac[1], &cc[1], &p[1], &ap[1]);
844 pap = xdot_(nx, ny, nz, &p[1], &ap[1]);
845 alpha = rhok2 / pap;
846 /* * */
847 /* * *** save rhok2 for next iteration *** */
848 rhok1 = rhok2;
849 /* * */
850 /* * *** update solution in direction p of length alpha *** */
851 xaxpy_(nx, ny, nz, &alpha, &p[1], &x[1]);
852 /* * */
853 /* * *** update residual *** */
854 d__1 = -alpha;
855 xaxpy_(nx, ny, nz, &d__1, &ap[1], &r__[1]);
856 /* * */
857 /* * *** some bookkeeping *** */
858 ++(*iters);
859 /* * */
860 /* * *** compute/check the current stopping test *** */
861 orsnrm = rsnrm;
862 if (*istop == 0) {
863 rsnrm = xnrm1_(nx, ny, nz, &r__[1]);
864 } else if (*istop == 1) {
865 rsnrm = xnrm1_(nx, ny, nz, &r__[1]);
866 } else if (*istop == 2) {
867 xcopy_(nx, ny, nz, &tru[1], &tmp[1]);
868 xaxpy_(nx, ny, nz, &c_b73, &x[1], &tmp[1]);
869 rsnrm = xnrm1_(nx, ny, nz, &tmp[1]);
870 } else if (*istop == 3) {
871 xcopy_(nx, ny, nz, &tru[1], &tmp[1]);
872 xaxpy_(nx, ny, nz, &c_b73, &x[1], &tmp[1]);
873 rsnrm = xnrm2_(nx, ny, nz, &tmp[1]);
874 } else if (*istop == 4) {
875 xcopy_(nx, ny, nz, &tru[1], &tmp[1]);
876 xaxpy_(nx, ny, nz, &c_b73, &x[1], &tmp[1]);
877 rsnrm = xnrm2_(nx, ny, nz, &tmp[1]);
878 } else if (*istop == 5) {
879 xcopy_(nx, ny, nz, &tru[1], &tmp[1]);
880 xaxpy_(nx, ny, nz, &c_b73, &x[1], &tmp[1]);
881 matvec_(nx, ny, nz, &ipc[1], &rpc[1], &ac[1], &cc[1], &tmp[1], &zk[1])
882 ;
883 rsnrm = sqrt(xdot_(nx, ny, nz, &tmp[1], &zk[1]));
884 } else {
885 s_wsle(&io___80);
886 do_lio(&c__9, &c__1, "% CGHSGO: bad istop value... ", (ftnlen)29);
887 e_wsle();
888 }
889 prtstp_(iok, iters, &rsnrm, &rsden, &orsnrm);
890 if (rsnrm / rsden <= *errtol) {
891 goto L99;
892 }
893 if (*iters >= *itmax) {
894 goto L99;
895 }
896 goto L30;
897 /* * */
898 /* * *** return and end *** */
899 L99:
900 return 0;
901 } /* cghsgo_ */
902
903