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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 8206 byte(s)
Setting up web subdirectory in repository
1 aw0a 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