1 /* linalg/condest.c
2 *
3 * Copyright (C) 2016 Patrick Alken
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include <config.h>
21 #include <stdlib.h>
22 #include <string.h>
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_vector.h>
25 #include <gsl/gsl_matrix.h>
26 #include <gsl/gsl_blas.h>
27 #include <gsl/gsl_linalg.h>
28
29 /*
30 * This module contains routines for estimating the condition number
31 * of matrices in the 1-norm. The algorithm is based on the paper,
32 *
33 * [1] N. J. Higham, "FORTRAN codes for estimating the one-norm of
34 * a real or complex matrix, with applications to condition estimation",
35 * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
36 */
37
38 static double condest_tri_norm1(CBLAS_UPLO_t Uplo, const gsl_matrix * A);
39 static int condest_tri_rcond(CBLAS_UPLO_t Uplo, const gsl_matrix * A,
40 double * rcond, gsl_vector * work);
41 static int condest_same_sign(const gsl_vector * x, const gsl_vector * y);
42 static int condest_invtriu(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);
43 static int condest_invtril(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);
44
45 int
gsl_linalg_tri_rcond(CBLAS_UPLO_t Uplo,const gsl_matrix * A,double * rcond,gsl_vector * work)46 gsl_linalg_tri_rcond(CBLAS_UPLO_t Uplo, const gsl_matrix * A, double * rcond, gsl_vector * work)
47 {
48 return condest_tri_rcond(Uplo, A, rcond, work);
49 }
50
51 /*
52 gsl_linalg_invnorm1()
53 Estimate the 1-norm of ||A^{-1}||, where A is a square
54 N-by-N matrix
55
56 Inputs: N - size of matrix
57 Ainvx - pointer to function which calculates:
58 x := A^{-1} x or x := A^{-t} x
59 params - parameters to pass to Ainvx
60 Ainvnorm - (output) estimate of ||A^{-1}||_1
61 work - workspace, length 3*N
62 */
63
64 int
gsl_linalg_invnorm1(const size_t N,int (* Ainvx)(CBLAS_TRANSPOSE_t TransA,gsl_vector * x,void * params),void * params,double * Ainvnorm,gsl_vector * work)65 gsl_linalg_invnorm1(const size_t N,
66 int (* Ainvx)(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params),
67 void * params, double * Ainvnorm, gsl_vector * work)
68 {
69 if (work->size != 3 * N)
70 {
71 GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
72 }
73 else
74 {
75 const size_t maxit = 5;
76 gsl_vector_view x = gsl_vector_subvector(work, 0, N);
77 gsl_vector_view v = gsl_vector_subvector(work, N, N);
78 gsl_vector_view xi = gsl_vector_subvector(work, 2*N, N);
79 double gamma, gamma_old, temp;
80 size_t i, k;
81
82 for (i = 0; i < N; ++i)
83 gsl_vector_set(&x.vector, i, 1.0 / (double) N);
84
85 /* compute v = A^{-1} x */
86 gsl_vector_memcpy(&v.vector, &x.vector);
87 (*Ainvx)(CblasNoTrans, &v.vector, params);
88
89 /* gamma = ||v||_1 */
90 gamma = gsl_blas_dasum(&v.vector);
91
92 /* xi = sign(v) */
93 for (i = 0; i < N; ++i)
94 {
95 double vi = gsl_vector_get(&v.vector, i);
96 gsl_vector_set(&xi.vector, i, GSL_SIGN(vi));
97 }
98
99 /* x = A^{-t} xi */
100 gsl_vector_memcpy(&x.vector, &xi.vector);
101 (*Ainvx)(CblasTrans, &x.vector, params);
102
103 for (k = 0; k < maxit; ++k)
104 {
105 size_t j = (size_t) gsl_blas_idamax(&x.vector);
106
107 /* v := A^{-1} e_j */
108 gsl_vector_set_zero(&v.vector);
109 gsl_vector_set(&v.vector, j, 1.0);
110 (*Ainvx)(CblasNoTrans, &v.vector, params);
111
112 gamma_old = gamma;
113 gamma = gsl_blas_dasum(&v.vector);
114
115 /* check for repeated sign vector (algorithm has converged) */
116 if (condest_same_sign(&v.vector, &xi.vector) || (gamma < gamma_old))
117 break;
118
119 /* xi = sign(v) */
120 for (i = 0; i < N; ++i)
121 {
122 double vi = gsl_vector_get(&v.vector, i);
123 gsl_vector_set(&xi.vector, i, GSL_SIGN(vi));
124 }
125
126 /* x = A^{-t} sign(v) */
127 gsl_vector_memcpy(&x.vector, &xi.vector);
128 (*Ainvx)(CblasTrans, &x.vector, params);
129 }
130
131 temp = 1.0; /* (-1)^i */
132 for (i = 0; i < N; ++i)
133 {
134 double term = 1.0 + (double) i / (N - 1.0);
135 gsl_vector_set(&x.vector, i, temp * term);
136 temp = -temp;
137 }
138
139 /* x := A^{-1} x */
140 (*Ainvx)(CblasNoTrans, &x.vector, params);
141
142 temp = 2.0 * gsl_blas_dasum(&x.vector) / (3.0 * N);
143 if (temp > gamma)
144 {
145 gsl_vector_memcpy(&v.vector, &x.vector);
146 gamma = temp;
147 }
148
149 *Ainvnorm = gamma;
150
151 return GSL_SUCCESS;
152 }
153 }
154
155 static int
condest_tri_rcond(CBLAS_UPLO_t Uplo,const gsl_matrix * A,double * rcond,gsl_vector * work)156 condest_tri_rcond(CBLAS_UPLO_t Uplo, const gsl_matrix * A, double * rcond, gsl_vector * work)
157 {
158 const size_t M = A->size1;
159 const size_t N = A->size2;
160
161 if (M != N)
162 {
163 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
164 }
165 else if (work->size != 3 * N)
166 {
167 GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
168 }
169 else
170 {
171 int status;
172 double Anorm = condest_tri_norm1(Uplo, A); /* ||A||_1 */
173 double Ainvnorm; /* ||A^{-1}||_1 */
174
175 *rcond = 0.0;
176
177 /* don't continue if matrix is singular */
178 if (Anorm == 0.0)
179 return GSL_SUCCESS;
180
181 /* estimate ||A^{-1}||_1 */
182 if (Uplo == CblasUpper)
183 status = gsl_linalg_invnorm1(N, condest_invtriu, (void *) A, &Ainvnorm, work);
184 else
185 status = gsl_linalg_invnorm1(N, condest_invtril, (void *) A, &Ainvnorm, work);
186
187 if (status)
188 return status;
189
190 if (Ainvnorm != 0.0)
191 *rcond = (1.0 / Anorm) / Ainvnorm;
192
193 return GSL_SUCCESS;
194 }
195 }
196
197 /* calculate 1 norm of triangular matrix */
198 static double
condest_tri_norm1(CBLAS_UPLO_t Uplo,const gsl_matrix * A)199 condest_tri_norm1(CBLAS_UPLO_t Uplo, const gsl_matrix * A)
200 {
201 const size_t N = A->size2;
202 double max = 0.0;
203 size_t i, j;
204
205 if (Uplo == CblasUpper)
206 {
207 for (j = 0; j < N; ++j)
208 {
209 double sum = 0.0;
210 for (i = 0; i <= j; ++i)
211 {
212 double Aij = gsl_matrix_get(A, i, j);
213 sum += fabs(Aij);
214 }
215
216 max = GSL_MAX(max, sum);
217 }
218 }
219 else
220 {
221 for (j = 0; j < N; ++j)
222 {
223 double sum = 0.0;
224 for (i = j; i < N; ++i)
225 {
226 double Aij = gsl_matrix_get(A, i, j);
227 sum += fabs(Aij);
228 }
229
230 max = GSL_MAX(max, sum);
231 }
232 }
233
234 return max;
235 }
236
237 /* return 1 if sign(x) = sign(y), 0 otherwise */
238 static int
condest_same_sign(const gsl_vector * x,const gsl_vector * y)239 condest_same_sign(const gsl_vector * x, const gsl_vector * y)
240 {
241 const size_t n = x->size;
242 size_t i;
243
244 for (i = 0; i < n; ++i)
245 {
246 double xi = gsl_vector_get(x, i);
247 double yi = gsl_vector_get(y, i);
248 if (GSL_SIGN(xi) != GSL_SIGN(yi))
249 return 0;
250 }
251
252 return 1;
253 }
254
255 /* x := A^{-1} x, A upper triangular */
256 static int
condest_invtriu(CBLAS_TRANSPOSE_t TransA,gsl_vector * x,void * params)257 condest_invtriu(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
258 {
259 gsl_matrix * A = (gsl_matrix *) params;
260 return gsl_blas_dtrsv(CblasUpper, TransA, CblasNonUnit, A, x);
261 }
262
263 /* x := A^{-1} x, A lower triangular */
264 static int
condest_invtril(CBLAS_TRANSPOSE_t TransA,gsl_vector * x,void * params)265 condest_invtril(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
266 {
267 gsl_matrix * A = (gsl_matrix *) params;
268 return gsl_blas_dtrsv(CblasLower, TransA, CblasNonUnit, A, x);
269 }
270
271 #ifndef GSL_DISABLE_DEPRECATED
272
273 int
gsl_linalg_tri_upper_rcond(const gsl_matrix * A,double * rcond,gsl_vector * work)274 gsl_linalg_tri_upper_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work)
275 {
276 int status = condest_tri_rcond(CblasUpper, A, rcond, work);
277 return status;
278 }
279
280 int
gsl_linalg_tri_lower_rcond(const gsl_matrix * A,double * rcond,gsl_vector * work)281 gsl_linalg_tri_lower_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work)
282 {
283 int status = condest_tri_rcond(CblasLower, A, rcond, work);
284 return status;
285 }
286
287 #endif
288