template_lapack_tptri.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_TPTRI_HEADER
00036 #define TEMPLATE_LAPACK_TPTRI_HEADER
00037 
00038 #include "template_lapack_common.h"
00039 
00040 template<class Treal>
00041 int template_lapack_tptri(const char *uplo, const char *diag, const integer *n, Treal *
00042         ap, 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        September 30, 1994   
00048 
00049 
00050     Purpose   
00051     =======   
00052 
00053     DTPTRI computes the inverse of a real upper or lower triangular   
00054     matrix A stored in packed format.   
00055 
00056     Arguments   
00057     =========   
00058 
00059     UPLO    (input) CHARACTER*1   
00060             = 'U':  A is upper triangular;   
00061             = 'L':  A is lower triangular.   
00062 
00063     DIAG    (input) CHARACTER*1   
00064             = 'N':  A is non-unit triangular;   
00065             = 'U':  A is unit triangular.   
00066 
00067     N       (input) INTEGER   
00068             The order of the matrix A.  N >= 0.   
00069 
00070     AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)   
00071             On entry, the upper or lower triangular matrix A, stored   
00072             columnwise in a linear array.  The j-th column of A is stored   
00073             in the array AP as follows:   
00074             if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;   
00075             if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n.   
00076             See below for further details.   
00077             On exit, the (triangular) inverse of the original matrix, in   
00078             the same packed storage format.   
00079 
00080     INFO    (output) INTEGER   
00081             = 0:  successful exit   
00082             < 0:  if INFO = -i, the i-th argument had an illegal value   
00083             > 0:  if INFO = i, A(i,i) is exactly zero.  The triangular   
00084                   matrix is singular and its inverse can not be computed.   
00085 
00086     Further Details   
00087     ===============   
00088 
00089     A triangular matrix A can be transferred to packed storage using one   
00090     of the following program segments:   
00091 
00092     UPLO = 'U':                      UPLO = 'L':   
00093 
00094           JC = 1                           JC = 1   
00095           DO 2 J = 1, N                    DO 2 J = 1, N   
00096              DO 1 I = 1, J                    DO 1 I = J, N   
00097                 AP(JC+I-1) = A(I,J)              AP(JC+I-J) = A(I,J)   
00098         1    CONTINUE                    1    CONTINUE   
00099              JC = JC + J                      JC = JC + N - J + 1   
00100         2 CONTINUE                       2 CONTINUE   
00101 
00102     =====================================================================   
00103 
00104 
00105        Test the input parameters.   
00106 
00107        Parameter adjustments */
00108     /* Table of constant values */
00109      integer c__1 = 1;
00110     
00111     /* System generated locals */
00112     integer i__1, i__2;
00113     /* Local variables */
00114      integer j;
00115      logical upper;
00116      integer jc, jj;
00117      integer jclast;
00118      logical nounit;
00119      Treal ajj;
00120 
00121 
00122     --ap;
00123 
00124     /* Initialization added by Elias to get rid of compiler warnings. */
00125     jclast = 0;
00126     /* Function Body */
00127     *info = 0;
00128     upper = template_blas_lsame(uplo, "U");
00129     nounit = template_blas_lsame(diag, "N");
00130     if (! upper && ! template_blas_lsame(uplo, "L")) {
00131         *info = -1;
00132     } else if (! nounit && ! template_blas_lsame(diag, "U")) {
00133         *info = -2;
00134     } else if (*n < 0) {
00135         *info = -3;
00136     }
00137     if (*info != 0) {
00138         i__1 = -(*info);
00139         template_blas_erbla("TPTRI ", &i__1);
00140         return 0;
00141     }
00142 
00143 /*     Check for singularity if non-unit. */
00144 
00145     if (nounit) {
00146         if (upper) {
00147             jj = 0;
00148             i__1 = *n;
00149             for (*info = 1; *info <= i__1; ++(*info)) {
00150                 jj += *info;
00151                 if (ap[jj] == 0.) {
00152                     return 0;
00153                 }
00154 /* L10: */
00155             }
00156         } else {
00157             jj = 1;
00158             i__1 = *n;
00159             for (*info = 1; *info <= i__1; ++(*info)) {
00160                 if (ap[jj] == 0.) {
00161                     return 0;
00162                 }
00163                 jj = jj + *n - *info + 1;
00164 /* L20: */
00165             }
00166         }
00167         *info = 0;
00168     }
00169 
00170     if (upper) {
00171 
00172 /*        Compute inverse of upper triangular matrix. */
00173 
00174         jc = 1;
00175         i__1 = *n;
00176         for (j = 1; j <= i__1; ++j) {
00177             if (nounit) {
00178                 ap[jc + j - 1] = 1. / ap[jc + j - 1];
00179                 ajj = -ap[jc + j - 1];
00180             } else {
00181                 ajj = -1.;
00182             }
00183 
00184 /*           Compute elements 1:j-1 of j-th column. */
00185 
00186             i__2 = j - 1;
00187             template_blas_tpmv("Upper", "No transpose", diag, &i__2, &ap[1], &ap[jc], &
00188                     c__1);
00189             i__2 = j - 1;
00190             template_blas_scal(&i__2, &ajj, &ap[jc], &c__1);
00191             jc += j;
00192 /* L30: */
00193         }
00194 
00195     } else {
00196 
00197 /*        Compute inverse of lower triangular matrix. */
00198 
00199         jc = *n * (*n + 1) / 2;
00200         for (j = *n; j >= 1; --j) {
00201             if (nounit) {
00202                 ap[jc] = 1. / ap[jc];
00203                 ajj = -ap[jc];
00204             } else {
00205                 ajj = -1.;
00206             }
00207             if (j < *n) {
00208 
00209 /*              Compute elements j+1:n of j-th column. */
00210 
00211                 i__1 = *n - j;
00212                 template_blas_tpmv("Lower", "No transpose", diag, &i__1, &ap[jclast], &ap[
00213                         jc + 1], &c__1);
00214                 i__1 = *n - j;
00215                 template_blas_scal(&i__1, &ajj, &ap[jc + 1], &c__1);
00216             }
00217             jclast = jc;
00218             jc = jc - *n + j - 2;
00219 /* L40: */
00220         }
00221     }
00222 
00223     return 0;
00224 
00225 /*     End of DTPTRI */
00226 
00227 } /* dtptri_ */
00228 
00229 #endif

Generated on Wed Nov 21 09:32:38 2012 for ergo by  doxygen 1.4.7