template_lapack_trti2.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_TRTI2_HEADER
00036 #define TEMPLATE_LAPACK_TRTI2_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_trti2(const char *uplo, const char *diag, const integer *n, Treal *
00041         a, const integer *lda, integer *info)
00042 {
00043 /*  -- LAPACK routine (version 3.0) --   
00044        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00045        Courant Institute, Argonne National Lab, and Rice University   
00046        February 29, 1992   
00047 
00048 
00049     Purpose   
00050     =======   
00051 
00052     DTRTI2 computes the inverse of a real upper or lower triangular   
00053     matrix.   
00054 
00055     This is the Level 2 BLAS version of the algorithm.   
00056 
00057     Arguments   
00058     =========   
00059 
00060     UPLO    (input) CHARACTER*1   
00061             Specifies whether the matrix A is upper or lower triangular.   
00062             = 'U':  Upper triangular   
00063             = 'L':  Lower triangular   
00064 
00065     DIAG    (input) CHARACTER*1   
00066             Specifies whether or not the matrix A is unit triangular.   
00067             = 'N':  Non-unit triangular   
00068             = 'U':  Unit triangular   
00069 
00070     N       (input) INTEGER   
00071             The order of the matrix A.  N >= 0.   
00072 
00073     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)   
00074             On entry, the triangular matrix A.  If UPLO = 'U', the   
00075             leading n by n upper triangular part of the array A contains   
00076             the upper triangular matrix, and the strictly lower   
00077             triangular part of A is not referenced.  If UPLO = 'L', the   
00078             leading n by n lower triangular part of the array A contains   
00079             the lower triangular matrix, and the strictly upper   
00080             triangular part of A is not referenced.  If DIAG = 'U', the   
00081             diagonal elements of A are also not referenced and are   
00082             assumed to be 1.   
00083 
00084             On exit, the (triangular) inverse of the original matrix, in   
00085             the same storage format.   
00086 
00087     LDA     (input) INTEGER   
00088             The leading dimension of the array A.  LDA >= max(1,N).   
00089 
00090     INFO    (output) INTEGER   
00091             = 0: successful exit   
00092             < 0: if INFO = -k, the k-th argument had an illegal value   
00093 
00094     =====================================================================   
00095 
00096 
00097        Test the input parameters.   
00098 
00099        Parameter adjustments */
00100     /* Table of constant values */
00101      integer c__1 = 1;
00102     
00103     /* System generated locals */
00104     integer a_dim1, a_offset, i__1, i__2;
00105     /* Local variables */
00106      integer j;
00107      logical upper;
00108      logical nounit;
00109      Treal ajj;
00110 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00111 
00112 
00113     a_dim1 = *lda;
00114     a_offset = 1 + a_dim1 * 1;
00115     a -= a_offset;
00116 
00117     /* Function Body */
00118     *info = 0;
00119     upper = template_blas_lsame(uplo, "U");
00120     nounit = template_blas_lsame(diag, "N");
00121     if (! upper && ! template_blas_lsame(uplo, "L")) {
00122         *info = -1;
00123     } else if (! nounit && ! template_blas_lsame(diag, "U")) {
00124         *info = -2;
00125     } else if (*n < 0) {
00126         *info = -3;
00127     } else if (*lda < maxMACRO(1,*n)) {
00128         *info = -5;
00129     }
00130     if (*info != 0) {
00131         i__1 = -(*info);
00132         template_blas_erbla("TRTI2 ", &i__1);
00133         return 0;
00134     }
00135 
00136     if (upper) {
00137 
00138 /*        Compute inverse of upper triangular matrix. */
00139 
00140         i__1 = *n;
00141         for (j = 1; j <= i__1; ++j) {
00142             if (nounit) {
00143                 a_ref(j, j) = 1. / a_ref(j, j);
00144                 ajj = -a_ref(j, j);
00145             } else {
00146                 ajj = -1.;
00147             }
00148 
00149 /*           Compute elements 1:j-1 of j-th column. */
00150 
00151             i__2 = j - 1;
00152             template_blas_trmv("Upper", "No transpose", diag, &i__2, &a[a_offset], lda, &
00153                     a_ref(1, j), &c__1);
00154             i__2 = j - 1;
00155             template_blas_scal(&i__2, &ajj, &a_ref(1, j), &c__1);
00156 /* L10: */
00157         }
00158     } else {
00159 
00160 /*        Compute inverse of lower triangular matrix. */
00161 
00162         for (j = *n; j >= 1; --j) {
00163             if (nounit) {
00164                 a_ref(j, j) = 1. / a_ref(j, j);
00165                 ajj = -a_ref(j, j);
00166             } else {
00167                 ajj = -1.;
00168             }
00169             if (j < *n) {
00170 
00171 /*              Compute elements j+1:n of j-th column. */
00172 
00173                 i__1 = *n - j;
00174                 template_blas_trmv("Lower", "No transpose", diag, &i__1, &a_ref(j + 1, j 
00175                         + 1), lda, &a_ref(j + 1, j), &c__1);
00176                 i__1 = *n - j;
00177                 template_blas_scal(&i__1, &ajj, &a_ref(j + 1, j), &c__1);
00178             }
00179 /* L20: */
00180         }
00181     }
00182 
00183     return 0;
00184 
00185 /*     End of DTRTI2 */
00186 
00187 } /* dtrti2_ */
00188 
00189 #undef a_ref
00190 
00191 
00192 #endif

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