ergo
template_lapack_pocon.h
Go to the documentation of this file.
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