1 #include "alberta_util.h"
2 #include "alberta_util_intern.h"
3
4 /*--------------------------------------------------------------------------*/
5 /*--- information output for linear solvers ------------------------------*/
6 /*--------------------------------------------------------------------------*/
7
start_info(const char * funcName,OEM_DATA * oem)8 void start_info(const char *funcName, OEM_DATA *oem)
9 {
10 oem->info = oem->info > 10 ? 10 : oem->info;
11
12 INFO(oem->info,6,"with tolerance %le", oem->tolerance);
13 if (oem->restart > 0)
14 PRINT_INFO(oem->info,6," and restart %d\n", oem->restart);
15 else
16 PRINT_INFO(oem->info,6,"\n");
17 INFO(oem->info,2,"iter. | residual | red.\n");
18 fflush(stdout);
19 return;
20 }
21
break_info(const char * funcName,OEM_DATA * oem,const char * reason,int iter,REAL res,REAL * ores,WORKSPACE * ws)22 void break_info(const char *funcName, OEM_DATA *oem, const char *reason,
23 int iter, REAL res, REAL *ores, WORKSPACE *ws)
24 {
25 if (*ores && *ores > 0)
26 INFO(oem->info,2,"%5d | %12.5le | %8.2le\n", iter, res, res/(*ores));
27 else
28 INFO(oem->info,2,"%5d | %12.5le |\n", iter);
29 INFO(oem->info,2,"stop due to: %s\n", reason);
30 fflush(stdout);
31
32 free_oem_workspace(ws, oem);
33 oem->residual = res;
34 }
35
solve_info(const char * funcName,OEM_DATA * oem,int iter,REAL res,REAL * ores,WORKSPACE * ws)36 int solve_info(const char *funcName, OEM_DATA *oem, int iter, REAL res,
37 REAL *ores, WORKSPACE *ws)
38 {
39 static int step[11] = {0, 1000, 500, 200, 100, 50, 20, 10, 5, 2, 1};
40
41 if (res <= oem->tolerance || (oem->info && (iter%step[oem->info] == 0))
42 || iter == oem->max_iter)
43 {
44 if (*ores)
45 {
46 if (*ores > 0.0)
47 {
48 REAL red = res/(*ores);
49 INFO(oem->info,2,"%5d | %12.5le | %8.2le\n", iter, res, red);
50 }
51 else
52 {
53 INFO(oem->info,2,"%5d | %12.5le | --------\n", iter, res);
54 }
55 *ores = res;
56 }
57 else
58 {
59 INFO(oem->info,2,"%5d | %12.5le |\n", iter, res);
60 }
61 }
62 oem->residual = res;
63
64 if (iter >= oem->max_iter || res <= oem->tolerance) {
65 if (res > oem->tolerance) {
66 INFO(oem->info,1,"tolerance %le not reached after %d iterations\n",
67 oem->tolerance, iter);
68 } else {
69 INFO(oem->info,6,"finished successfully with %d iterations\n",iter);
70 }
71 fflush(stdout);
72 free_oem_workspace(ws, oem);
73 return 1;
74 }
75
76 fflush(stdout);
77 return 0;
78 }
79
80 /*--------------------------------------------------------------------------*/
81 /*--- checking of workspace, reallocation of workspace if neccesary ------*/
82 /*--------------------------------------------------------------------------*/
83
check_workspace(const char * funcName,const char * file,int line,size_t size,WORKSPACE * ws)84 WORKSPACE *check_workspace(const char *funcName, const char *file, int line,
85 size_t size, WORKSPACE *ws)
86 {
87 if (!ws)
88 {
89 ws = GET_WORKSPACE(size*sizeof(REAL));
90 }
91 else if (size*sizeof(REAL) > ws->size)
92 {
93 WARNING("need workspace for %d REALs\n", size);
94 WARNING("reallocating workspace of length %d\n", size*sizeof(REAL));
95 REALLOC_WORKSPACE(ws, size*sizeof(REAL));
96 }
97 return(ws);
98 }
99
free_oem_workspace(WORKSPACE * ws,OEM_DATA * oem)100 void free_oem_workspace(WORKSPACE *ws, OEM_DATA *oem)
101 {
102 if (ws != oem->ws)
103 FREE_WORKSPACE(ws);
104 return;
105 }
106