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