ascend4/solver/linutils.c

/*
 *  lin_utils: Ascend Linear Algebra Utilities
 *  by Kirk Andre' Abbott
 *  Created: 12 March 1995
 *  Version: $Revision: 1.5 $
 *  Version control file: $RCSfile: linutils.c,v $
 *  Date last modified: $Date: 2000/01/25 02:27:01 $
 *  Last modified by: $Author: ballan $
 *
 *  This file is part of the SLV solver.
 *
 *  Copyright (C) 1995 Kirk Andre' Abbott
 *
 *  The SLV solver is free software; you can redistribute
 *  it and/or modify it under the terms of the GNU General Public License as
 *  published by the Free Software Foundation; either version 2 of the
 *  License, or (at your option) any later version.
 *
 *  The SLV solver is distributed in hope that it will be
 *  useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License along with
 *  the program; if not, write to the Free Software Foundation, Inc., 675
 *  Mass Ave, Cambridge, MA 02139 USA.  Check the file named COPYING.
 *  COPYING is found in ../compiler.
 */

#include <stdio.h>
#include <math.h>
#include "utilities/ascConfig.h"
#include "compiler/compiler.h"
#include "utilities/ascMalloc.h"
#include "general/tm_time.h"
#include "solver/mtx.h"
#include "solver/linsol.h"
#include "solver/linsolqr.h"
#include "solver/linutils.h"

/*
 * The A1 norm is:
 *   ||A||\1 = max (i<=j<=n) SUM\1..m |Aij|
 * Because we have the mtx_sum_sqrs_in_col operator, we will instead
 * do:
 *   ||A||\1 = max (i<=j<=n) SUM\1..m (Aij)^2.0
 *
 */

double linutils_A_1_norm(mtx_matrix_t mtx, mtx_region_t *reg)

{
  mtx_region_t mtxreg;
  double value, maxvalue = 0;
  int col;

  if (mtx==NULL)
    return 1.0e20;

  if (reg==NULL) {
    mtxreg.row.low = mtxreg.col.low = 0;
    mtxreg.row.high = mtxreg.col.high = mtx_order(mtx);
  }
  else{
    mtxreg = *reg;
  }
  for (col = mtxreg.col.low; col < mtxreg.col.high; col++) {
    value = mtx_sum_sqrs_in_col(mtx, col, &mtxreg.row);
    maxvalue = MAX(maxvalue,value);
  }
  return maxvalue;
}


/*
 * The A\infinity norm is:
 *   ||A||\infinity = max (i<=i<=m) SUM\1..n |Aij|
 * where i,m refers to rows and j,n refers to columns.
 * Because we have the mtx_sum_sqrs_in_row operator, we will instead
 * do:
 *   ||A||\infinity =   max      SUM   (Aij)^2.0
 *           (i<=i<=m) (1..n)
 */

double linutils_A_infinity_norm(mtx_matrix_t mtx, mtx_region_t *reg)

{
  mtx_region_t mtxreg;
  double value, maxvalue = 0;
  int row;

  if (mtx==NULL)
    return 1.0e20;

  if (reg==NULL) {
    mtxreg.row.low = mtxreg.col.low = 0;
    mtxreg.row.high = mtxreg.col.high = mtx_order(mtx);
  }
  else{
    mtxreg = *reg;
  }
  for (row = mtxreg.row.low; row < mtxreg.row.high; row++) {
    value = mtx_sum_sqrs_in_row(mtx, row, &mtxreg.col);
    maxvalue = MAX(maxvalue,value);
  }
  return maxvalue;
}

/*
 * The Frobenius norm is given as:
 *
 * ||A||\F =   ( sum     sum   Aij^2 )^ 1/2
 *            1<=j<=n  1<=i<=m
 */
double linutils_A_Frobenius_norm(mtx_matrix_t mtx, mtx_region_t *reg)
{
  mtx_region_t mtxreg;
  double value;
  /* double maxvalue = 0; */
  int row;

  if (mtx==NULL)
    return 1.0e20;

  if (reg==NULL) {
    mtxreg.row.low = mtxreg.col.low = 0;
    mtxreg.row.high = mtxreg.col.high = mtx_order(mtx);
  }
  else{
    mtxreg = *reg;
  }

  value = 0.0;
  for (row = mtxreg.row.low; row < mtxreg.row.high; row++) {
    value += mtx_sum_sqrs_in_row(mtx, row, &mtxreg.col);
  }
  return sqrt(value);
}

/*
 * This function has this funny name, because I am using
 * a very inelegant method to do this computation, but I
 * am using the linsolqr system primitives. Ben Allan has
 * much more efficient ways of doing these computations,
 * in fact, he gets them for free at the end of his QR
 * analyses. It assumes that the linear system provided
 * has *already* been factored. We then compute the inverse
 * explicitly, compute the Frobenius norms for both the
 * inverse and the original matrix and return:
 *               -1
 * || A || . || A ||
 *       F         F
 * Warning !! this function is neither fast nor memory
 * efficient. !!
 */
double linutils_A_condqr_kaa(linsolqr_system_t lin_sys,
                 mtx_matrix_t mtx,
                 mtx_region_t *reg)
{
  mtx_matrix_t factors;
  mtx_region_t mtxreg;
  double *solution = NULL, *rhs = NULL;
  double a_norm, a_inverse_norm = 0.0;
  int rank,capacity;
  int j,k;
  double time;

  if (mtx==NULL)
    return 1.0e20;

  if (reg==NULL) {
    mtxreg.row.low = mtxreg.col.low = 0;
    mtxreg.row.high = mtxreg.col.high = mtx_symbolic_rank(mtx);
  }
  else{
    mtxreg = *reg;
  }

  /*
   * Let us do the expensive part first, which is to calculate,
   * A inverse column by column. While are doing this we will
   * compute and accumulated the 2 norm for each column.
   */
  factors = linsolqr_get_matrix(lin_sys);
  capacity = mtx_capacity(factors);
  solution = (double *)asccalloc(capacity,sizeof(double));
  rhs = (double *)asccalloc(capacity,sizeof(double));
  linsolqr_add_rhs(lin_sys,rhs,FALSE);

  rank = mtx_symbolic_rank(factors);
  time = tm_cpu_time();
  for (j=0;j<rank;j++) {
    rhs[j] = 1.0;
    linsolqr_rhs_was_changed(lin_sys,rhs);
    linsolqr_solve(lin_sys,rhs);
    linsolqr_copy_solution(lin_sys,rhs,solution);

    for (k=0;k<capacity;k++) {/* dense vector 2 norm */
      a_inverse_norm += solution[k]*solution[k];
      solution[k] = 0.0;
    }
    rhs[j] = 0.0;
  }
  time = tm_cpu_time() - time;
  FPRINTF(stderr,"Time to compute explicit inverse by linsolqr = %g\n",time);
  a_inverse_norm = sqrt(a_inverse_norm);
  linsolqr_remove_rhs(lin_sys,rhs);

  /*
   * We should now have a_inverse_norm
   * All we need to do now is to calculate a_norm and
   * return the result: a_norm * a_inverse_norm.
   */

  a_norm = linutils_A_Frobenius_norm(mtx, &mtxreg);
  FPRINTF(stderr,"A norm = %g\tA_inverse norm = %g; condition # = %g\n",
      a_norm,a_inverse_norm,a_norm*a_inverse_norm);

  if (solution) ascfree(solution);
  if (rhs) ascfree(rhs);

  return (a_norm*a_inverse_norm);
}


/*
 * Same function as above, but using linsol primitives
 * rather than linsolqr primitives. This function is transient
 * and is here for temporary comparison only, of the time it
 * takes to do the inverse computations.
 */
double linutils_A_cond_kaa(linsol_system_t lin_sys,
               mtx_matrix_t mtx,
               mtx_region_t *reg)
{
  mtx_matrix_t factors;
  mtx_region_t mtxreg;
  double *solution = NULL, *rhs = NULL;
  double a_norm, a_inverse_norm = 0.0;
  int rank,capacity;
  int j,k;
  double time;

  if (mtx==NULL)
    return 1.0e20;

  if (reg==NULL) {
    mtxreg.row.low = mtxreg.col.low = 0;
    mtxreg.row.high = mtxreg.col.high = mtx_symbolic_rank(mtx);
  }
  else{
    mtxreg = *reg;
  }

  /*
   * Let us do the expensive part first, which is to calculate,
   * A inverse column by column. While are doing this we will
   * compute and accumulated the 2 norm for each column.
   */
  factors = linsol_get_matrix(lin_sys);
  capacity = mtx_capacity(factors);
  solution = (double *)asccalloc(capacity,sizeof(double));
  rhs = (double *)asccalloc(capacity,sizeof(double));
  linsol_add_rhs(lin_sys,rhs,FALSE);

  rank = mtx_symbolic_rank(factors);
  time = tm_cpu_time();
  for (j=0;j<rank;j++) {
    rhs[j] = 1.0;
    linsol_rhs_was_changed(lin_sys,rhs);
    linsol_solve(lin_sys,rhs);
    linsol_copy_solution(lin_sys,rhs,solution);

    for (k=0;k<capacity;k++) {/* dense vector 2 norm */
      a_inverse_norm += solution[k]*solution[k];
      solution[k] = 0.0;
    }
    rhs[j] = 0.0;
  }
  time = tm_cpu_time() - time;
  FPRINTF(stderr,"Time to compute explicit inverse by linsol = %g\n",time);
  a_inverse_norm = sqrt(a_inverse_norm);
  linsol_remove_rhs(lin_sys,rhs);

  /*
   * We should now have a_inverse_norm
   * All we need to do now is to calculate a_norm and
   * return the result: a_norm * a_inverse_norm.
   */

  a_norm = linutils_A_Frobenius_norm(mtx, &mtxreg);
  FPRINTF(stderr,"A norm = %g\tA_inverse norm = %g; condition # = %g\n",
      a_norm,a_inverse_norm,a_norm*a_inverse_norm);

  if (solution) ascfree(solution);
  if (rhs) ascfree(rhs);

  return (a_norm*a_inverse_norm);
}

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