ergo
template_lapack_pocon.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /* This file belongs to the template_lapack part of the Ergo source
31 * code. The source files in the template_lapack directory are modified
32 * versions of files originally distributed as CLAPACK, see the
33 * Copyright/license notice in the file template_lapack/COPYING.
34 */
35
36
37#ifndef TEMPLATE_LAPACK_POCON_HEADER
38#define TEMPLATE_LAPACK_POCON_HEADER
39
40
41template<class Treal>
42int template_lapack_pocon(const char *uplo, const integer *n, const Treal *a, const integer *
43 lda, const Treal *anorm, Treal *rcond, Treal *work, integer *
44 iwork, integer *info)
45{
46/* -- LAPACK routine (version 3.0) --
47 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48 Courant Institute, Argonne National Lab, and Rice University
49 March 31, 1993
50
51
52 Purpose
53 =======
54
55 DPOCON estimates the reciprocal of the condition number (in the
56 1-norm) of a real symmetric positive definite matrix using the
57 Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF.
58
59 An estimate is obtained for norm(inv(A)), and the reciprocal of the
60 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
61
62 Arguments
63 =========
64
65 UPLO (input) CHARACTER*1
66 = 'U': Upper triangle of A is stored;
67 = 'L': Lower triangle of A is stored.
68
69 N (input) INTEGER
70 The order of the matrix A. N >= 0.
71
72 A (input) DOUBLE PRECISION array, dimension (LDA,N)
73 The triangular factor U or L from the Cholesky factorization
74 A = U**T*U or A = L*L**T, as computed by DPOTRF.
75
76 LDA (input) INTEGER
77 The leading dimension of the array A. LDA >= max(1,N).
78
79 ANORM (input) DOUBLE PRECISION
80 The 1-norm (or infinity-norm) of the symmetric matrix A.
81
82 RCOND (output) DOUBLE PRECISION
83 The reciprocal of the condition number of the matrix A,
84 computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
85 estimate of the 1-norm of inv(A) computed in this routine.
86
87 WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
88
89 IWORK (workspace) INTEGER array, dimension (N)
90
91 INFO (output) INTEGER
92 = 0: successful exit
93 < 0: if INFO = -i, the i-th argument had an illegal value
94
95 =====================================================================
96
97
98 Test the input parameters.
99
100 Parameter adjustments */
101 /* Table of constant values */
102 integer c__1 = 1;
103
104 /* System generated locals */
105 integer a_dim1, a_offset, i__1;
106 Treal d__1;
107 /* Local variables */
108 integer kase;
109 Treal scale;
110 logical upper;
111 integer ix;
112 Treal scalel;
113 Treal scaleu;
114 Treal ainvnm;
115 char normin[1];
116 Treal smlnum;
117
118
119 a_dim1 = *lda;
120 a_offset = 1 + a_dim1 * 1;
121 a -= a_offset;
122 --work;
123 --iwork;
124
125 /* Function Body */
126 *info = 0;
127 upper = template_blas_lsame(uplo, "U");
128 if (! upper && ! template_blas_lsame(uplo, "L")) {
129 *info = -1;
130 } else if (*n < 0) {
131 *info = -2;
132 } else if (*lda < maxMACRO(1,*n)) {
133 *info = -4;
134 } else if (*anorm < 0.) {
135 *info = -5;
136 }
137 if (*info != 0) {
138 i__1 = -(*info);
139 template_blas_erbla("POCON ", &i__1);
140 return 0;
141 }
142
143/* Quick return if possible */
144
145 *rcond = 0.;
146 if (*n == 0) {
147 *rcond = 1.;
148 return 0;
149 } else if (*anorm == 0.) {
150 return 0;
151 }
152
153 smlnum = template_lapack_lamch("Safe minimum", (Treal)0);
154
155/* Estimate the 1-norm of inv(A). */
156
157 kase = 0;
158 *(unsigned char *)normin = 'N';
159L10:
160 template_lapack_lacon(n, &work[*n + 1], &work[1], &iwork[1], &ainvnm, &kase);
161 if (kase != 0) {
162 if (upper) {
163
164/* Multiply by inv(U'). */
165
166 template_lapack_latrs("Upper", "Transpose", "Non-unit", normin, n, &a[a_offset],
167 lda, &work[1], &scalel, &work[(*n << 1) + 1], info);
168 *(unsigned char *)normin = 'Y';
169
170/* Multiply by inv(U). */
171
172 template_lapack_latrs("Upper", "No transpose", "Non-unit", normin, n, &a[
173 a_offset], lda, &work[1], &scaleu, &work[(*n << 1) + 1],
174 info);
175 } else {
176
177/* Multiply by inv(L). */
178
179 template_lapack_latrs("Lower", "No transpose", "Non-unit", normin, n, &a[
180 a_offset], lda, &work[1], &scalel, &work[(*n << 1) + 1],
181 info);
182 *(unsigned char *)normin = 'Y';
183
184/* Multiply by inv(L'). */
185
186 template_lapack_latrs("Lower", "Transpose", "Non-unit", normin, n, &a[a_offset],
187 lda, &work[1], &scaleu, &work[(*n << 1) + 1], info);
188 }
189
190/* Multiply by 1/SCALE if doing so will not cause overflow. */
191
192 scale = scalel * scaleu;
193 if (scale != 1.) {
194 ix = template_blas_idamax(n, &work[1], &c__1);
195 if (scale < (d__1 = work[ix], absMACRO(d__1)) * smlnum || scale == 0.)
196 {
197 goto L20;
198 }
199 template_lapack_rscl(n, &scale, &work[1], &c__1);
200 }
201 goto L10;
202 }
203
204/* Compute the estimate of the reciprocal condition number. */
205
206 if (ainvnm != 0.) {
207 *rcond = 1. / ainvnm / *anorm;
208 }
209
210L20:
211 return 0;
212
213/* End of DPOCON */
214
215} /* dpocon_ */
216
217#endif
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
int integer
Definition: template_blas_common.h:40
#define absMACRO(x)
Definition: template_blas_common.h:47
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
integer template_blas_idamax(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_idamax.h:42
int template_lapack_lacon(const integer *n, Treal *v, Treal *x, integer *isgn, Treal *est, integer *kase)
Definition: template_lapack_lacon.h:42
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
int template_lapack_latrs(const char *uplo, const char *trans, const char *diag, const char *normin, const integer *n, const Treal *a, const integer *lda, Treal *x, Treal *scale, Treal *cnorm, integer *info)
Definition: template_lapack_latrs.h:42
int template_lapack_pocon(const char *uplo, const integer *n, const Treal *a, const integer *lda, const Treal *anorm, Treal *rcond, Treal *work, integer *iwork, integer *info)
Definition: template_lapack_pocon.h:42
int template_lapack_rscl(const integer *n, const Treal *sa, Treal *sx, const integer *incx)
Definition: template_lapack_rscl.h:42