ergo
template_lapack_pocon.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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 
41 template<class Treal>
42 int 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';
159 L10:
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 
210 L20:
211  return 0;
212 
213 /* End of DPOCON */
214 
215 } /* dpocon_ */
216 
217 #endif
#define absMACRO(x)
Definition: template_blas_common.h:47
integer template_blas_idamax(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_idamax.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 integer
Definition: template_blas_common.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
int template_lapack_lacon(const integer *n, Treal *v, Treal *x, integer *isgn, Treal *est, integer *kase)
Definition: template_lapack_lacon.h:42
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
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
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
int template_lapack_rscl(const integer *n, const Treal *sa, Treal *sx, const integer *incx)
Definition: template_lapack_rscl.h:42
bool logical
Definition: template_blas_common.h:41
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46