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