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	/*
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