1 /*
2 * lu_condest.c
3 *
4 * Copyright (C) 2016-2018 ERGO-Code
5 *
6 * LINPACK condition number estimate
7 *
8 */
9
10 #include "lu_internal.h"
11
12 /*
13 * lu_condest()
14 *
15 * Given m-by-m matrix U such that U[perm,perm] is upper triangular,
16 * return estimate for 1-norm condition number of U.
17 * If @norm is not NULL, it holds the 1-norm of the matrix on return.
18 * If @norminv is not NULL, it holds the estimated 1-norm of the inverse on
19 * return.
20 *
21 * The other function arguments are the same as in lu_normest().
22 *
23 */
lu_condest(lu_int m,const lu_int * Ubegin,const lu_int * Ui,const double * Ux,const double * pivot,const lu_int * perm,int upper,double * work,double * norm,double * norminv)24 double lu_condest(
25 lu_int m, const lu_int *Ubegin, const lu_int *Ui, const double *Ux,
26 const double *pivot, const lu_int *perm, int upper, double *work,
27 double *norm, double *norminv)
28 {
29 lu_int j, p;
30 double Unorm, Uinvnorm;
31
32 /* compute 1-norm of U */
33 Unorm = 0;
34 for (j = 0; j < m; j++)
35 {
36 double colsum = pivot ? fabs(pivot[j]) : 1;
37 for (p = Ubegin[j]; Ui[p] >= 0; p++)
38 colsum += fabs(Ux[p]);
39 Unorm = fmax(Unorm, colsum);
40 }
41
42 /* estimate 1-norm of U^{-1} */
43 Uinvnorm = lu_normest(m, Ubegin, Ui, Ux, pivot, perm, upper, work);
44
45 if (norm) *norm = Unorm;
46 if (norminv) *norminv = Uinvnorm;
47
48 return Unorm * Uinvnorm;
49 }
50
51 /*
52 * lu_normest()
53 *
54 * Given m-by-m matrix U such that U[perm,perm] is triangular,
55 * estimate 1-norm of U^{-1} by computing
56 *
57 * U'x = b, Uy = x, normest = max{norm(y)_1/norm(x)_1, norm(x)_inf},
58 *
59 * where the entries of b are +/-1 chosen dynamically to make x large.
60 * The method is described in [1].
61 *
62 * @Ubegin, @Ui, @Ux matrix U in compressed column format without pivots,
63 * columns are terminated by a negative index
64 * @pivot pivot elements by column index of U; NULL if unit pivots
65 * @perm permutation to triangular form; NULL if identity
66 * @upper nonzero if permuted matrix is upper triangular; zero if lower
67 * @work size m workspace, uninitialized on entry/return
68 *
69 * Return: estimate for 1-norm of U^{-1}
70 *
71 * [1] I. Duff, A. Erisman, J. Reid, "Direct Methods for Sparse Matrices"
72 *
73 */
lu_normest(lu_int m,const lu_int * Ubegin,const lu_int * Ui,const double * Ux,const double * pivot,const lu_int * perm,int upper,double * work)74 double lu_normest(
75 lu_int m, const lu_int *Ubegin, const lu_int *Ui, const double *Ux,
76 const double *pivot, const lu_int *perm, int upper, double *work)
77 {
78 lu_int i, j, k, kbeg, kend, kinc, p;
79 double x1norm, xinfnorm, y1norm, temp;
80
81 x1norm = 0;
82 xinfnorm = 0;
83 if (upper)
84 {
85 kbeg = 0; kend = m; kinc = 1;
86 }
87 else
88 {
89 kbeg = m-1; kend = -1; kinc = -1;
90 }
91 for (k = kbeg; k != kend; k += kinc)
92 {
93 j = perm ? perm[k] : k;
94 temp = 0;
95 for (p = Ubegin[j]; (i = Ui[p]) >= 0; p++)
96 temp -= work[i] * Ux[p];
97 temp += temp >= 0 ? 1 : -1; /* choose b[i] = 1 or b[i] = -1 */
98 if (pivot) temp /= pivot[j];
99 work[j] = temp;
100 x1norm += fabs(temp);
101 xinfnorm = fmax(xinfnorm, fabs(temp));
102 }
103
104 y1norm = 0;
105 if (upper)
106 {
107 kbeg = m-1; kend = -1; kinc = -1;
108 }
109 else
110 {
111 kbeg = 0; kend = m; kinc = 1;
112 }
113 for (k = kbeg; k != kend; k += kinc)
114 {
115 j = perm ? perm[k] : k;
116 if (pivot) work[j] /= pivot[j];
117 temp = work[j];
118 for (p = Ubegin[j]; (i = Ui[p]) >= 0; p++)
119 work[i] -= temp * Ux[p];
120 y1norm += fabs(temp);
121 }
122
123 return fmax(y1norm/x1norm, xinfnorm);
124 }
125