1 /* -----------------------------------------------------------------
2 * Programmer(s): Allan Taylor, Alan Hindmarsh and
3 * Radu Serban @ LLNL
4 * -----------------------------------------------------------------
5 * SUNDIALS Copyright Start
6 * Copyright (c) 2002-2021, Lawrence Livermore National Security
7 * and Southern Methodist University.
8 * All rights reserved.
9 *
10 * See the top-level LICENSE and NOTICE files for details.
11 *
12 * SPDX-License-Identifier: BSD-3-Clause
13 * SUNDIALS Copyright End
14 * -----------------------------------------------------------------
15 *
16 * This example loops through the available iterative linear solvers:
17 * SPGMR, SPBCG and SPTFQMR.
18 *
19 * Example problem for IDA: 2D heat equation, serial, GMRES.
20 *
21 * This example solves a discretized 2D heat equation problem.
22 * This version loops through the Krylov solvers Spgmr, Spbcg
23 * and Sptfqmr.
24 *
25 * The DAE system solved is a spatial discretization of the PDE
26 * du/dt = d^2u/dx^2 + d^2u/dy^2
27 * on the unit square. The boundary condition is u = 0 on all edges.
28 * Initial conditions are given by u = 16 x (1 - x) y (1 - y). The
29 * PDE is treated with central differences on a uniform M x M grid.
30 * The values of u at the interior points satisfy ODEs, and
31 * equations u = 0 at the boundaries are appended, to form a DAE
32 * system of size N = M^2. Here M = 10.
33 *
34 * The system is solved with IDA using the following Krylov
35 * linear solvers: SPGMR, SPBCG and SPTFQMR. The
36 * preconditioner uses the diagonal elements of the Jacobian only.
37 * Routines for preconditioning, required by SP*, are supplied
38 * here. The constraints u >= 0 are posed for all components. Output
39 * is taken at t = 0, .01, .02, .04,..., 10.24.
40 * -----------------------------------------------------------------*/
41
42 #include <stdio.h>
43 #include <stdlib.h>
44 #include <math.h>
45
46 #include <idas/idas.h> /* main integrator header file */
47 #include <sunlinsol/sunlinsol_spgmr.h> /* access to SPGMR SUNLinearSolver */
48 #include <sunlinsol/sunlinsol_spbcgs.h> /* access to SPBCGS SUNLinearSolver */
49 #include <sunlinsol/sunlinsol_sptfqmr.h> /* access to SPTFQMR SUNLinearSolver */
50 #include <nvector/nvector_serial.h> /* serial N_Vector types, fct. and macros */
51 #include <sundials/sundials_types.h> /* definition of realtype */
52
53 /* helpful macros */
54
55 #ifndef SQRT
56 #if defined(SUNDIALS_DOUBLE_PRECISION)
57 #define SQRT(x) (sqrt((x)))
58 #elif defined(SUNDIALS_SINGLE_PRECISION)
59 #define SQRT(x) (sqrtf((x)))
60 #elif defined(SUNDIALS_EXTENDED_PRECISION)
61 #define SQRT(x) (sqrtl((x)))
62 #endif
63 #endif
64
65 /* Problem Constants */
66
67 #define NOUT 11
68 #define MGRID 10
69 #define NEQ MGRID*MGRID
70 #define ZERO RCONST(0.0)
71 #define ONE RCONST(1.0)
72 #define TWO RCONST(2.0)
73 #define FOUR RCONST(4.0)
74
75 /* Linear Solver Loop Constants */
76
77 #define USE_SPGMR 0
78 #define USE_SPBCG 1
79 #define USE_SPTFQMR 2
80
81 /* User data type */
82
83 typedef struct {
84 sunindextype mm; /* number of grid points */
85 realtype dx;
86 realtype coeff;
87 N_Vector pp; /* vector of prec. diag. elements */
88 } *UserData;
89
90 /* Prototypes for functions called by IDA */
91
92 int resHeat(realtype tres, N_Vector uu, N_Vector up,
93 N_Vector resval, void *user_data);
94
95 int PsetupHeat(realtype tt,
96 N_Vector uu, N_Vector up, N_Vector rr,
97 realtype c_j, void *user_data);
98
99 int PsolveHeat(realtype tt,
100 N_Vector uu, N_Vector up, N_Vector rr,
101 N_Vector rvec, N_Vector zvec,
102 realtype c_j, realtype delta, void *user_data);
103
104 /* Prototypes for private functions */
105
106 static int SetInitialProfile(UserData data, N_Vector uu, N_Vector up,
107 N_Vector res);
108 static void PrintHeader(realtype rtol, realtype atol, int linsolver);
109 static void PrintOutput(void *mem, realtype t, N_Vector uu, int linsolver);
110 static int check_retval(void *returnvalue, const char *funcname, int opt);
111
112 /*
113 *--------------------------------------------------------------------
114 * MAIN PROGRAM
115 *--------------------------------------------------------------------
116 */
117
main(int argc,char * argv[])118 int main(int argc, char* argv[])
119 {
120 void *mem;
121 UserData data;
122 N_Vector uu, up, constraints, res;
123 int retval, iout, linsolver;
124 realtype rtol, atol, t0, t1, tout, tret;
125 long int netf, ncfn, ncfl;
126 SUNLinearSolver LS;
127 int nrmfactor; /* LS norm conversion factor flag */
128 realtype nrmfac; /* LS norm conversion factor */
129
130 mem = NULL;
131 data = NULL;
132 uu = NULL;
133 up = NULL;
134 constraints = NULL;
135 res = NULL;
136 LS = NULL;
137 nrmfactor = 0;
138
139 /* Retrieve the command-line options */
140 if (argc > 1) nrmfactor = atoi(argv[1]);
141
142 /* Allocate N-vectors and the user data structure. */
143
144 uu = N_VNew_Serial(NEQ);
145 if(check_retval((void *)uu, "N_VNew_Serial", 0)) return(1);
146
147 up = N_VNew_Serial(NEQ);
148 if(check_retval((void *)up, "N_VNew_Serial", 0)) return(1);
149
150 res = N_VNew_Serial(NEQ);
151 if(check_retval((void *)res, "N_VNew_Serial", 0)) return(1);
152
153 constraints = N_VNew_Serial(NEQ);
154 if(check_retval((void *)constraints, "N_VNew_Serial", 0)) return(1);
155
156 data = (UserData) malloc(sizeof *data);
157 data->pp = NULL;
158 if(check_retval((void *)data, "malloc", 2)) return(1);
159
160 /* Assign parameters in the user data structure. */
161
162 data->mm = MGRID;
163 data->dx = ONE/(MGRID-ONE);
164 data->coeff = ONE/(data->dx * data->dx);
165 data->pp = N_VNew_Serial(NEQ);
166 if(check_retval((void *)data->pp, "N_VNew_Serial", 0)) return(1);
167
168 /* Initialize uu, up. */
169
170 SetInitialProfile(data, uu, up, res);
171
172 /* Set constraints to all 1's for nonnegative solution values. */
173
174 N_VConst(ONE, constraints);
175
176 /* Assign various parameters. */
177
178 t0 = ZERO;
179 t1 = RCONST(0.01);
180 rtol = ZERO;
181 atol = RCONST(1.0e-3);
182
183 /* Call IDACreate and IDAMalloc to initialize solution */
184
185 mem = IDACreate();
186 if(check_retval((void *)mem, "IDACreate", 0)) return(1);
187
188 retval = IDASetUserData(mem, data);
189 if(check_retval(&retval, "IDASetUserData", 1)) return(1);
190
191 retval = IDASetConstraints(mem, constraints);
192 if(check_retval(&retval, "IDASetConstraints", 1)) return(1);
193 N_VDestroy(constraints);
194
195 retval = IDAInit(mem, resHeat, t0, uu, up);
196 if(check_retval(&retval, "IDAInit", 1)) return(1);
197
198 retval = IDASStolerances(mem, rtol, atol);
199 if(check_retval(&retval, "IDASStolerances", 1)) return(1);
200
201 /* START: Loop through SPGMR, SPBCG and SPTFQMR linear solver modules */
202 for (linsolver = 0; linsolver < 3; ++linsolver) {
203
204 if (linsolver != 0) {
205
206 /* Re-initialize uu, up. */
207 SetInitialProfile(data, uu, up, res);
208
209 /* Re-initialize IDA */
210 retval = IDAReInit(mem, t0, uu, up);
211 if (check_retval(&retval, "IDAReInit", 1)) return(1);
212
213 }
214
215 /* Free previous linear solver and attach a new linear solver module */
216 SUNLinSolFree(LS);
217
218 switch(linsolver) {
219
220 /* (a) SPGMR */
221 case(USE_SPGMR):
222
223 /* Print header */
224 printf(" -------");
225 printf(" \n| SPGMR |\n");
226 printf(" -------\n");
227
228 /* Call SUNLinSol_SPGMR to specify the linear solver SPGMR with
229 left preconditioning and the default maximum Krylov dimension */
230 LS = SUNLinSol_SPGMR(uu, PREC_LEFT, 0);
231 if(check_retval((void *)LS, "SUNLinSol_SPGMR", 0)) return(1);
232
233 /* Attach the linear solver */
234 retval = IDASetLinearSolver(mem, LS, NULL);
235 if(check_retval(&retval, "IDASetLinearSolver", 1)) return 1;
236
237 break;
238
239 /* (b) SPBCG */
240 case(USE_SPBCG):
241
242 /* Print header */
243 printf(" -------");
244 printf(" \n| SPBCGS |\n");
245 printf(" -------\n");
246
247 /* Call SUNLinSol_SPBCGS to specify the linear solver SPBCGS with
248 left preconditioning and the default maximum Krylov dimension */
249 LS = SUNLinSol_SPBCGS(uu, PREC_LEFT, 0);
250 if(check_retval((void *)LS, "SUNLinSol_SPBCGS", 0)) return(1);
251
252 /* Attach the linear solver */
253 retval = IDASetLinearSolver(mem, LS, NULL);
254 if(check_retval(&retval, "IDASetLinearSolver", 1)) return 1;
255
256 break;
257
258 /* (c) SPTFQMR */
259 case(USE_SPTFQMR):
260
261 /* Print header */
262 printf(" ---------");
263 printf(" \n| SPTFQMR |\n");
264 printf(" ---------\n");
265
266 /* Call SUNLinSol_SPTFQMR to specify the linear solver SPTFQMR with
267 left preconditioning and the default maximum Krylov dimension */
268 LS = SUNLinSol_SPTFQMR(uu, PREC_LEFT, 0);
269 if(check_retval((void *)LS, "SUNLinSol_SPTFQMR", 0)) return(1);
270
271 /* Attach the linear solver */
272 retval = IDASetLinearSolver(mem, LS, NULL);
273 if(check_retval(&retval, "IDASetLinearSolver", 1)) return 1;
274
275 break;
276
277 }
278
279 /* Specify preconditioner */
280 retval = IDASetPreconditioner(mem, PsetupHeat, PsolveHeat);
281 if(check_retval(&retval, "IDASetPreconditioner", 1)) return(1);
282
283 /* Set the linear solver tolerance conversion factor */
284 switch(nrmfactor) {
285
286 case(1):
287 /* use the square root of the vector length */
288 nrmfac = SQRT((realtype)NEQ);
289 break;
290 case(2):
291 /* compute with dot product */
292 nrmfac = -ONE;
293 break;
294 default:
295 /* use the default */
296 nrmfac = ZERO;
297 break;
298 }
299
300 retval = IDASetLSNormFactor(mem, nrmfac);
301 if (check_retval(&retval, "IDASetLSNormFactor", 1)) return(1);
302
303 /* Print output heading. */
304 PrintHeader(rtol, atol, linsolver);
305
306 /* Print output table heading, and initial line of table. */
307
308 printf("\n Output Summary (umax = max-norm of solution) \n\n");
309 printf(" time umax k nst nni nje nre nreLS h npe nps\n" );
310 printf("----------------------------------------------------------------------\n");
311
312 /* Loop over output times, call IDASolve, and print results. */
313
314 for (tout = t1,iout = 1; iout <= NOUT ; iout++, tout *= TWO) {
315 retval = IDASolve(mem, tout, &tret, uu, up, IDA_NORMAL);
316 if(check_retval(&retval, "IDASolve", 1)) return(1);
317 PrintOutput(mem, tret, uu, linsolver);
318 }
319
320 /* Print remaining counters. */
321 retval = IDAGetNumErrTestFails(mem, &netf);
322 check_retval(&retval, "IDAGetNumErrTestFails", 1);
323
324 retval = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
325 check_retval(&retval, "IDAGetNumNonlinSolvConvFails", 1);
326
327 retval = IDAGetNumLinConvFails(mem, &ncfl);
328 check_retval(&retval, "IDAGetNumLinConvFails", 1);
329
330 printf("\nError test failures = %ld\n", netf);
331 printf("Nonlinear convergence failures = %ld\n", ncfn);
332 printf("Linear convergence failures = %ld\n", ncfl);
333
334 if (linsolver < 2)
335 printf("\n======================================================================\n\n");
336
337 } /* END: Loop through SPGMR, SPBCG and SPTFQMR linear solver modules */
338
339 /* Free Memory */
340
341 IDAFree(&mem);
342 SUNLinSolFree(LS);
343
344 N_VDestroy(uu);
345 N_VDestroy(up);
346 N_VDestroy(res);
347
348 N_VDestroy(data->pp);
349 free(data);
350
351 return(0);
352 }
353
354 /*
355 *--------------------------------------------------------------------
356 * FUNCTIONS CALLED BY IDA
357 *--------------------------------------------------------------------
358 */
359
360 /*
361 * resHeat: heat equation system residual function (user-supplied)
362 * This uses 5-point central differencing on the interior points, and
363 * includes algebraic equations for the boundary values.
364 * So for each interior point, the residual component has the form
365 * res_i = u'_i - (central difference)_i
366 * while for each boundary point, it is res_i = u_i.
367 */
368
resHeat(realtype tt,N_Vector uu,N_Vector up,N_Vector rr,void * user_data)369 int resHeat(realtype tt,
370 N_Vector uu, N_Vector up, N_Vector rr,
371 void *user_data)
372 {
373 sunindextype i, j, offset, loc, mm;
374 realtype *uu_data, *up_data, *rr_data, coeff, dif1, dif2;
375 UserData data;
376
377 uu_data = N_VGetArrayPointer(uu);
378 up_data = N_VGetArrayPointer(up);
379 rr_data = N_VGetArrayPointer(rr);
380
381 data = (UserData) user_data;
382
383 coeff = data->coeff;
384 mm = data->mm;
385
386 /* Initialize rr to uu, to take care of boundary equations. */
387 N_VScale(ONE, uu, rr);
388
389 /* Loop over interior points; set res = up - (central difference). */
390 for (j = 1; j < MGRID-1; j++) {
391 offset = mm*j;
392 for (i = 1; i < mm-1; i++) {
393 loc = offset + i;
394 dif1 = uu_data[loc-1] + uu_data[loc+1] - TWO * uu_data[loc];
395 dif2 = uu_data[loc-mm] + uu_data[loc+mm] - TWO * uu_data[loc];
396 rr_data[loc]= up_data[loc] - coeff * ( dif1 + dif2 );
397 }
398 }
399
400 return(0);
401 }
402
403 /*
404 * PsetupHeat: setup for diagonal preconditioner.
405 *
406 * The optional user-supplied functions PsetupHeat and
407 * PsolveHeat together must define the left preconditoner
408 * matrix P approximating the system Jacobian matrix
409 * J = dF/du + cj*dF/du'
410 * (where the DAE system is F(t,u,u') = 0), and solve the linear
411 * systems P z = r. This is done in this case by keeping only
412 * the diagonal elements of the J matrix above, storing them as
413 * inverses in a vector pp, when computed in PsetupHeat, for
414 * subsequent use in PsolveHeat.
415 *
416 * In this instance, only cj and data (user data structure, with
417 * pp etc.) are used from the PsetupdHeat argument list.
418 */
419
PsetupHeat(realtype tt,N_Vector uu,N_Vector up,N_Vector rr,realtype c_j,void * user_data)420 int PsetupHeat(realtype tt,
421 N_Vector uu, N_Vector up, N_Vector rr,
422 realtype c_j, void *user_data)
423 {
424
425 sunindextype i, j, offset, loc, mm;
426 realtype *ppv, pelinv;
427 UserData data;
428
429 data = (UserData) user_data;
430 ppv = N_VGetArrayPointer(data->pp);
431 mm = data->mm;
432
433 /* Initialize the entire vector to 1., then set the interior points to the
434 correct value for preconditioning. */
435 N_VConst(ONE,data->pp);
436
437 /* Compute the inverse of the preconditioner diagonal elements. */
438 pelinv = ONE/(c_j + FOUR*data->coeff);
439
440 for (j = 1; j < mm-1; j++) {
441 offset = mm * j;
442 for (i = 1; i < mm-1; i++) {
443 loc = offset + i;
444 ppv[loc] = pelinv;
445 }
446 }
447
448 return(0);
449 }
450
451 /*
452 * PsolveHeat: solve preconditioner linear system.
453 * This routine multiplies the input vector rvec by the vector pp
454 * containing the inverse diagonal Jacobian elements (previously
455 * computed in PrecondHeateq), returning the result in zvec.
456 */
457
PsolveHeat(realtype tt,N_Vector uu,N_Vector up,N_Vector rr,N_Vector rvec,N_Vector zvec,realtype c_j,realtype delta,void * user_data)458 int PsolveHeat(realtype tt,
459 N_Vector uu, N_Vector up, N_Vector rr,
460 N_Vector rvec, N_Vector zvec,
461 realtype c_j, realtype delta, void *user_data)
462 {
463 UserData data;
464 data = (UserData) user_data;
465 N_VProd(data->pp, rvec, zvec);
466 return(0);
467 }
468
469 /*
470 *--------------------------------------------------------------------
471 * PRIVATE FUNCTIONS
472 *--------------------------------------------------------------------
473 */
474
475 /*
476 * SetInitialProfile: routine to initialize u and up vectors.
477 */
478
SetInitialProfile(UserData data,N_Vector uu,N_Vector up,N_Vector res)479 static int SetInitialProfile(UserData data, N_Vector uu, N_Vector up,
480 N_Vector res)
481 {
482 sunindextype mm, mm1, i, j, offset, loc;
483 realtype xfact, yfact, *udata, *updata;
484
485 mm = data->mm;
486
487 udata = N_VGetArrayPointer(uu);
488 updata = N_VGetArrayPointer(up);
489
490 /* Initialize uu on all grid points. */
491 mm1 = mm - 1;
492 for (j = 0; j < mm; j++) {
493 yfact = data->dx * j;
494 offset = mm*j;
495 for (i = 0;i < mm; i++) {
496 xfact = data->dx * i;
497 loc = offset + i;
498 udata[loc] = RCONST(16.0) * xfact * (ONE - xfact) * yfact * (ONE - yfact);
499 }
500 }
501
502 /* Initialize up vector to 0. */
503 N_VConst(ZERO, up);
504
505 /* resHeat sets res to negative of ODE RHS values at interior points. */
506 resHeat(ZERO, uu, up, res, data);
507
508 /* Copy -res into up to get correct interior initial up values. */
509 N_VScale(-ONE, res, up);
510
511 /* Set up at boundary points to zero. */
512 for (j = 0; j < mm; j++) {
513 offset = mm*j;
514 for (i = 0; i < mm; i++) {
515 loc = offset + i;
516 if (j == 0 || j == mm1 || i == 0 || i == mm1 ) updata[loc] = ZERO;
517 }
518 }
519
520 return(0);
521 }
522
523 /*
524 * Print first lines of output (problem description)
525 */
526
PrintHeader(realtype rtol,realtype atol,int linsolver)527 static void PrintHeader(realtype rtol, realtype atol, int linsolver)
528 {
529 printf("\nidasKrylovDemo_ls: Heat equation, serial example problem for IDA\n");
530 printf(" Discretized heat equation on 2D unit square.\n");
531 printf(" Zero boundary conditions,");
532 printf(" polynomial initial conditions.\n");
533 printf(" Mesh dimensions: %d x %d", MGRID, MGRID);
534 printf(" Total system size: %d\n\n", NEQ);
535 #if defined(SUNDIALS_EXTENDED_PRECISION)
536 printf("Tolerance parameters: rtol = %Lg atol = %Lg\n", rtol, atol);
537 #elif defined(SUNDIALS_DOUBLE_PRECISION)
538 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
539 #else
540 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
541 #endif
542 printf("Constraints set to force all solution components >= 0. \n");
543
544 switch(linsolver) {
545
546 case(USE_SPGMR):
547 printf("Linear solver: SPGMR, preconditioner using diagonal elements. \n");
548 break;
549
550 case(USE_SPBCG):
551 printf("Linear solver: SPBCG, preconditioner using diagonal elements. \n");
552 break;
553
554 case(USE_SPTFQMR):
555 printf("Linear solver: SPTFQMR, preconditioner using diagonal elements. \n");
556 break;
557 }
558 }
559
560 /*
561 * PrintOutput: print max norm of solution and current solver statistics
562 */
563
PrintOutput(void * mem,realtype t,N_Vector uu,int linsolver)564 static void PrintOutput(void *mem, realtype t, N_Vector uu, int linsolver)
565 {
566 realtype hused, umax;
567 long int nst, nni, nje, nre, nreLS, nli, npe, nps;
568 int kused, retval;
569
570 umax = N_VMaxNorm(uu);
571
572 retval = IDAGetLastOrder(mem, &kused);
573 check_retval(&retval, "IDAGetLastOrder", 1);
574 retval = IDAGetNumSteps(mem, &nst);
575 check_retval(&retval, "IDAGetNumSteps", 1);
576 retval = IDAGetNumNonlinSolvIters(mem, &nni);
577 check_retval(&retval, "IDAGetNumNonlinSolvIters", 1);
578 retval = IDAGetNumResEvals(mem, &nre);
579 check_retval(&retval, "IDAGetNumResEvals", 1);
580 retval = IDAGetLastStep(mem, &hused);
581 check_retval(&retval, "IDAGetLastStep", 1);
582
583 retval = IDAGetNumJtimesEvals(mem, &nje);
584 check_retval(&retval, "IDAGetNumJtimesEvals", 1);
585 retval = IDAGetNumLinIters(mem, &nli);
586 check_retval(&retval, "IDAGetNumLinIters", 1);
587 retval = IDAGetNumLinResEvals(mem, &nreLS);
588 check_retval(&retval, "IDAGetNumLinResEvals", 1);
589 retval = IDAGetNumPrecEvals(mem, &npe);
590 check_retval(&retval, "IDAGetNumPrecEvals", 1);
591 retval = IDAGetNumPrecSolves(mem, &nps);
592 check_retval(&retval, "IDAGetNumPrecSolves", 1);
593
594 #if defined(SUNDIALS_EXTENDED_PRECISION)
595 printf(" %5.2Lf %13.5Le %d %3ld %3ld %3ld %4ld %4ld %9.2Le %3ld %3ld\n",
596 t, umax, kused, nst, nni, nje, nre, nreLS, hused, npe, nps);
597 #elif defined(SUNDIALS_DOUBLE_PRECISION)
598 printf(" %5.2f %13.5e %d %3ld %3ld %3ld %4ld %4ld %9.2e %3ld %3ld\n",
599 t, umax, kused, nst, nni, nje, nre, nreLS, hused, npe, nps);
600 #else
601 printf(" %5.2f %13.5e %d %3ld %3ld %3ld %4ld %4ld %9.2e %3ld %3ld\n",
602 t, umax, kused, nst, nni, nje, nre, nreLS, hused, npe, nps);
603 #endif
604 }
605
606 /*
607 * Check function return value...
608 * opt == 0 means SUNDIALS function allocates memory so check if
609 * returned NULL pointer
610 * opt == 1 means SUNDIALS function returns an integer value so check if
611 * retval < 0
612 * opt == 2 means function allocates memory so check if returned
613 * NULL pointer
614 */
615
check_retval(void * returnvalue,const char * funcname,int opt)616 static int check_retval(void *returnvalue, const char *funcname, int opt)
617 {
618 int *retval;
619
620 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
621 if (opt == 0 && returnvalue == NULL) {
622 fprintf(stderr,
623 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
624 funcname);
625 return(1);
626 } else if (opt == 1) {
627 /* Check if retval < 0 */
628 retval = (int *) returnvalue;
629 if (*retval < 0) {
630 fprintf(stderr,
631 "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
632 funcname, *retval);
633 return(1);
634 }
635 } else if (opt == 2 && returnvalue == NULL) {
636 /* Check if function returned NULL pointer - no memory allocated */
637 fprintf(stderr,
638 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
639 funcname);
640 return(1);
641 }
642
643 return(0);
644 }
645