/[ascend]/trunk/ascend4/solver/linutils.c
ViewVC logotype

Contents of /trunk/ascend4/solver/linutils.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 8206 byte(s)
Setting up web subdirectory in repository
1 /*
2 * lin_utils: Ascend Linear Algebra Utilities
3 * by Kirk Andre' Abbott
4 * Created: 12 March 1995
5 * Version: $Revision: 1.5 $
6 * Version control file: $RCSfile: linutils.c,v $
7 * Date last modified: $Date: 2000/01/25 02:27:01 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the SLV solver.
11 *
12 * Copyright (C) 1995 Kirk Andre' Abbott
13 *
14 * The SLV solver is free software; you can redistribute
15 * it and/or modify it under the terms of the GNU General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
18 *
19 * The SLV solver is distributed in hope that it will be
20 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22 * General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License along with
25 * the program; if not, write to the Free Software Foundation, Inc., 675
26 * Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
27 * COPYING is found in ../compiler.
28 */
29
30 #include <stdio.h>
31 #include <math.h>
32 #include "utilities/ascConfig.h"
33 #include "compiler/compiler.h"
34 #include "utilities/ascMalloc.h"
35 #include "general/tm_time.h"
36 #include "solver/mtx.h"
37 #include "solver/linsol.h"
38 #include "solver/linsolqr.h"
39 #include "solver/linutils.h"
40
41 /*
42 * The A1 norm is:
43 * ||A||\1 = max (i<=j<=n) SUM\1..m |Aij|
44 * Because we have the mtx_sum_sqrs_in_col operator, we will instead
45 * do:
46 * ||A||\1 = max (i<=j<=n) SUM\1..m (Aij)^2.0
47 *
48 */
49
50 double linutils_A_1_norm(mtx_matrix_t mtx, mtx_region_t *reg)
51
52 {
53 mtx_region_t mtxreg;
54 double value, maxvalue = 0;
55 int col;
56
57 if (mtx==NULL)
58 return 1.0e20;
59
60 if (reg==NULL) {
61 mtxreg.row.low = mtxreg.col.low = 0;
62 mtxreg.row.high = mtxreg.col.high = mtx_order(mtx);
63 }
64 else{
65 mtxreg = *reg;
66 }
67 for (col = mtxreg.col.low; col < mtxreg.col.high; col++) {
68 value = mtx_sum_sqrs_in_col(mtx, col, &mtxreg.row);
69 maxvalue = MAX(maxvalue,value);
70 }
71 return maxvalue;
72 }
73
74
75 /*
76 * The A\infinity norm is:
77 * ||A||\infinity = max (i<=i<=m) SUM\1..n |Aij|
78 * where i,m refers to rows and j,n refers to columns.
79 * Because we have the mtx_sum_sqrs_in_row operator, we will instead
80 * do:
81 * ||A||\infinity = max SUM (Aij)^2.0
82 * (i<=i<=m) (1..n)
83 */
84
85 double linutils_A_infinity_norm(mtx_matrix_t mtx, mtx_region_t *reg)
86
87 {
88 mtx_region_t mtxreg;
89 double value, maxvalue = 0;
90 int row;
91
92 if (mtx==NULL)
93 return 1.0e20;
94
95 if (reg==NULL) {
96 mtxreg.row.low = mtxreg.col.low = 0;
97 mtxreg.row.high = mtxreg.col.high = mtx_order(mtx);
98 }
99 else{
100 mtxreg = *reg;
101 }
102 for (row = mtxreg.row.low; row < mtxreg.row.high; row++) {
103 value = mtx_sum_sqrs_in_row(mtx, row, &mtxreg.col);
104 maxvalue = MAX(maxvalue,value);
105 }
106 return maxvalue;
107 }
108
109 /*
110 * The Frobenius norm is given as:
111 *
112 * ||A||\F = ( sum sum Aij^2 )^ 1/2
113 * 1<=j<=n 1<=i<=m
114 */
115 double linutils_A_Frobenius_norm(mtx_matrix_t mtx, mtx_region_t *reg)
116 {
117 mtx_region_t mtxreg;
118 double value;
119 /* double maxvalue = 0; */
120 int row;
121
122 if (mtx==NULL)
123 return 1.0e20;
124
125 if (reg==NULL) {
126 mtxreg.row.low = mtxreg.col.low = 0;
127 mtxreg.row.high = mtxreg.col.high = mtx_order(mtx);
128 }
129 else{
130 mtxreg = *reg;
131 }
132
133 value = 0.0;
134 for (row = mtxreg.row.low; row < mtxreg.row.high; row++) {
135 value += mtx_sum_sqrs_in_row(mtx, row, &mtxreg.col);
136 }
137 return sqrt(value);
138 }
139
140 /*
141 * This function has this funny name, because I am using
142 * a very inelegant method to do this computation, but I
143 * am using the linsolqr system primitives. Ben Allan has
144 * much more efficient ways of doing these computations,
145 * in fact, he gets them for free at the end of his QR
146 * analyses. It assumes that the linear system provided
147 * has *already* been factored. We then compute the inverse
148 * explicitly, compute the Frobenius norms for both the
149 * inverse and the original matrix and return:
150 * -1
151 * || A || . || A ||
152 * F F
153 * Warning !! this function is neither fast nor memory
154 * efficient. !!
155 */
156 double linutils_A_condqr_kaa(linsolqr_system_t lin_sys,
157 mtx_matrix_t mtx,
158 mtx_region_t *reg)
159 {
160 mtx_matrix_t factors;
161 mtx_region_t mtxreg;
162 double *solution = NULL, *rhs = NULL;
163 double a_norm, a_inverse_norm = 0.0;
164 int rank,capacity;
165 int j,k;
166 double time;
167
168 if (mtx==NULL)
169 return 1.0e20;
170
171 if (reg==NULL) {
172 mtxreg.row.low = mtxreg.col.low = 0;
173 mtxreg.row.high = mtxreg.col.high = mtx_symbolic_rank(mtx);
174 }
175 else{
176 mtxreg = *reg;
177 }
178
179 /*
180 * Let us do the expensive part first, which is to calculate,
181 * A inverse column by column. While are doing this we will
182 * compute and accumulated the 2 norm for each column.
183 */
184 factors = linsolqr_get_matrix(lin_sys);
185 capacity = mtx_capacity(factors);
186 solution = (double *)asccalloc(capacity,sizeof(double));
187 rhs = (double *)asccalloc(capacity,sizeof(double));
188 linsolqr_add_rhs(lin_sys,rhs,FALSE);
189
190 rank = mtx_symbolic_rank(factors);
191 time = tm_cpu_time();
192 for (j=0;j<rank;j++) {
193 rhs[j] = 1.0;
194 linsolqr_rhs_was_changed(lin_sys,rhs);
195 linsolqr_solve(lin_sys,rhs);
196 linsolqr_copy_solution(lin_sys,rhs,solution);
197
198 for (k=0;k<capacity;k++) {/* dense vector 2 norm */
199 a_inverse_norm += solution[k]*solution[k];
200 solution[k] = 0.0;
201 }
202 rhs[j] = 0.0;
203 }
204 time = tm_cpu_time() - time;
205 FPRINTF(stderr,"Time to compute explicit inverse by linsolqr = %g\n",time);
206 a_inverse_norm = sqrt(a_inverse_norm);
207 linsolqr_remove_rhs(lin_sys,rhs);
208
209 /*
210 * We should now have a_inverse_norm
211 * All we need to do now is to calculate a_norm and
212 * return the result: a_norm * a_inverse_norm.
213 */
214
215 a_norm = linutils_A_Frobenius_norm(mtx, &mtxreg);
216 FPRINTF(stderr,"A norm = %g\tA_inverse norm = %g; condition # = %g\n",
217 a_norm,a_inverse_norm,a_norm*a_inverse_norm);
218
219 if (solution) ascfree(solution);
220 if (rhs) ascfree(rhs);
221
222 return (a_norm*a_inverse_norm);
223 }
224
225
226 /*
227 * Same function as above, but using linsol primitives
228 * rather than linsolqr primitives. This function is transient
229 * and is here for temporary comparison only, of the time it
230 * takes to do the inverse computations.
231 */
232 double linutils_A_cond_kaa(linsol_system_t lin_sys,
233 mtx_matrix_t mtx,
234 mtx_region_t *reg)
235 {
236 mtx_matrix_t factors;
237 mtx_region_t mtxreg;
238 double *solution = NULL, *rhs = NULL;
239 double a_norm, a_inverse_norm = 0.0;
240 int rank,capacity;
241 int j,k;
242 double time;
243
244 if (mtx==NULL)
245 return 1.0e20;
246
247 if (reg==NULL) {
248 mtxreg.row.low = mtxreg.col.low = 0;
249 mtxreg.row.high = mtxreg.col.high = mtx_symbolic_rank(mtx);
250 }
251 else{
252 mtxreg = *reg;
253 }
254
255 /*
256 * Let us do the expensive part first, which is to calculate,
257 * A inverse column by column. While are doing this we will
258 * compute and accumulated the 2 norm for each column.
259 */
260 factors = linsol_get_matrix(lin_sys);
261 capacity = mtx_capacity(factors);
262 solution = (double *)asccalloc(capacity,sizeof(double));
263 rhs = (double *)asccalloc(capacity,sizeof(double));
264 linsol_add_rhs(lin_sys,rhs,FALSE);
265
266 rank = mtx_symbolic_rank(factors);
267 time = tm_cpu_time();
268 for (j=0;j<rank;j++) {
269 rhs[j] = 1.0;
270 linsol_rhs_was_changed(lin_sys,rhs);
271 linsol_solve(lin_sys,rhs);
272 linsol_copy_solution(lin_sys,rhs,solution);
273
274 for (k=0;k<capacity;k++) {/* dense vector 2 norm */
275 a_inverse_norm += solution[k]*solution[k];
276 solution[k] = 0.0;
277 }
278 rhs[j] = 0.0;
279 }
280 time = tm_cpu_time() - time;
281 FPRINTF(stderr,"Time to compute explicit inverse by linsol = %g\n",time);
282 a_inverse_norm = sqrt(a_inverse_norm);
283 linsol_remove_rhs(lin_sys,rhs);
284
285 /*
286 * We should now have a_inverse_norm
287 * All we need to do now is to calculate a_norm and
288 * return the result: a_norm * a_inverse_norm.
289 */
290
291 a_norm = linutils_A_Frobenius_norm(mtx, &mtxreg);
292 FPRINTF(stderr,"A norm = %g\tA_inverse norm = %g; condition # = %g\n",
293 a_norm,a_inverse_norm,a_norm*a_inverse_norm);
294
295 if (solution) ascfree(solution);
296 if (rhs) ascfree(rhs);
297
298 return (a_norm*a_inverse_norm);
299 }
300

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22