ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00028 /* This file belongs to the template_lapack part of the Ergo source 00029 * code. The source files in the template_lapack directory are modified 00030 * versions of files originally distributed as CLAPACK, see the 00031 * Copyright/license notice in the file template_lapack/COPYING. 00032 */ 00033 00034 00035 #ifndef TEMPLATE_LAPACK_POCON_HEADER 00036 #define TEMPLATE_LAPACK_POCON_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_pocon(const char *uplo, const integer *n, const Treal *a, const integer * 00041 lda, const Treal *anorm, Treal *rcond, Treal *work, integer * 00042 iwork, integer *info) 00043 { 00044 /* -- LAPACK routine (version 3.0) -- 00045 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00046 Courant Institute, Argonne National Lab, and Rice University 00047 March 31, 1993 00048 00049 00050 Purpose 00051 ======= 00052 00053 DPOCON estimates the reciprocal of the condition number (in the 00054 1-norm) of a real symmetric positive definite matrix using the 00055 Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF. 00056 00057 An estimate is obtained for norm(inv(A)), and the reciprocal of the 00058 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). 00059 00060 Arguments 00061 ========= 00062 00063 UPLO (input) CHARACTER*1 00064 = 'U': Upper triangle of A is stored; 00065 = 'L': Lower triangle of A is stored. 00066 00067 N (input) INTEGER 00068 The order of the matrix A. N >= 0. 00069 00070 A (input) DOUBLE PRECISION array, dimension (LDA,N) 00071 The triangular factor U or L from the Cholesky factorization 00072 A = U**T*U or A = L*L**T, as computed by DPOTRF. 00073 00074 LDA (input) INTEGER 00075 The leading dimension of the array A. LDA >= max(1,N). 00076 00077 ANORM (input) DOUBLE PRECISION 00078 The 1-norm (or infinity-norm) of the symmetric matrix A. 00079 00080 RCOND (output) DOUBLE PRECISION 00081 The reciprocal of the condition number of the matrix A, 00082 computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an 00083 estimate of the 1-norm of inv(A) computed in this routine. 00084 00085 WORK (workspace) DOUBLE PRECISION array, dimension (3*N) 00086 00087 IWORK (workspace) INTEGER array, dimension (N) 00088 00089 INFO (output) INTEGER 00090 = 0: successful exit 00091 < 0: if INFO = -i, the i-th argument had an illegal value 00092 00093 ===================================================================== 00094 00095 00096 Test the input parameters. 00097 00098 Parameter adjustments */ 00099 /* Table of constant values */ 00100 integer c__1 = 1; 00101 00102 /* System generated locals */ 00103 integer a_dim1, a_offset, i__1; 00104 Treal d__1; 00105 /* Local variables */ 00106 integer kase; 00107 Treal scale; 00108 logical upper; 00109 integer ix; 00110 Treal scalel; 00111 Treal scaleu; 00112 Treal ainvnm; 00113 char normin[1]; 00114 Treal smlnum; 00115 00116 00117 a_dim1 = *lda; 00118 a_offset = 1 + a_dim1 * 1; 00119 a -= a_offset; 00120 --work; 00121 --iwork; 00122 00123 /* Function Body */ 00124 *info = 0; 00125 upper = template_blas_lsame(uplo, "U"); 00126 if (! upper && ! template_blas_lsame(uplo, "L")) { 00127 *info = -1; 00128 } else if (*n < 0) { 00129 *info = -2; 00130 } else if (*lda < maxMACRO(1,*n)) { 00131 *info = -4; 00132 } else if (*anorm < 0.) { 00133 *info = -5; 00134 } 00135 if (*info != 0) { 00136 i__1 = -(*info); 00137 template_blas_erbla("POCON ", &i__1); 00138 return 0; 00139 } 00140 00141 /* Quick return if possible */ 00142 00143 *rcond = 0.; 00144 if (*n == 0) { 00145 *rcond = 1.; 00146 return 0; 00147 } else if (*anorm == 0.) { 00148 return 0; 00149 } 00150 00151 smlnum = template_lapack_lamch("Safe minimum", (Treal)0); 00152 00153 /* Estimate the 1-norm of inv(A). */ 00154 00155 kase = 0; 00156 *(unsigned char *)normin = 'N'; 00157 L10: 00158 template_lapack_lacon(n, &work[*n + 1], &work[1], &iwork[1], &ainvnm, &kase); 00159 if (kase != 0) { 00160 if (upper) { 00161 00162 /* Multiply by inv(U'). */ 00163 00164 template_lapack_latrs("Upper", "Transpose", "Non-unit", normin, n, &a[a_offset], 00165 lda, &work[1], &scalel, &work[(*n << 1) + 1], info); 00166 *(unsigned char *)normin = 'Y'; 00167 00168 /* Multiply by inv(U). */ 00169 00170 template_lapack_latrs("Upper", "No transpose", "Non-unit", normin, n, &a[ 00171 a_offset], lda, &work[1], &scaleu, &work[(*n << 1) + 1], 00172 info); 00173 } else { 00174 00175 /* Multiply by inv(L). */ 00176 00177 template_lapack_latrs("Lower", "No transpose", "Non-unit", normin, n, &a[ 00178 a_offset], lda, &work[1], &scalel, &work[(*n << 1) + 1], 00179 info); 00180 *(unsigned char *)normin = 'Y'; 00181 00182 /* Multiply by inv(L'). */ 00183 00184 template_lapack_latrs("Lower", "Transpose", "Non-unit", normin, n, &a[a_offset], 00185 lda, &work[1], &scaleu, &work[(*n << 1) + 1], info); 00186 } 00187 00188 /* Multiply by 1/SCALE if doing so will not cause overflow. */ 00189 00190 scale = scalel * scaleu; 00191 if (scale != 1.) { 00192 ix = template_blas_idamax(n, &work[1], &c__1); 00193 if (scale < (d__1 = work[ix], absMACRO(d__1)) * smlnum || scale == 0.) 00194 { 00195 goto L20; 00196 } 00197 template_lapack_rscl(n, &scale, &work[1], &c__1); 00198 } 00199 goto L10; 00200 } 00201 00202 /* Compute the estimate of the reciprocal condition number. */ 00203 00204 if (ainvnm != 0.) { 00205 *rcond = 1. / ainvnm / *anorm; 00206 } 00207 00208 L20: 00209 return 0; 00210 00211 /* End of DPOCON */ 00212 00213 } /* dpocon_ */ 00214 00215 #endif