template_lapack_trtri.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_TRTRI_HEADER
00036 #define TEMPLATE_LAPACK_TRTRI_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_trtri(const char *uplo, const char *diag,
00041                           const integer *n, Treal *a,
00042                           const integer *lda, 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     DTRTRI computes the inverse of a real upper or lower triangular   
00054     matrix A.   
00055 
00056     This is the Level 3 BLAS version of the algorithm.   
00057 
00058     Arguments   
00059     =========   
00060 
00061     UPLO    (input) CHARACTER*1   
00062             = 'U':  A is upper triangular;   
00063             = 'L':  A is lower triangular.   
00064 
00065     DIAG    (input) CHARACTER*1   
00066             = 'N':  A is non-unit triangular;   
00067             = 'U':  A is unit triangular.   
00068 
00069     N       (input) INTEGER   
00070             The order of the matrix A.  N >= 0.   
00071 
00072     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)   
00073             On entry, the triangular matrix A.  If UPLO = 'U', the   
00074             leading N-by-N upper triangular part of the array A contains   
00075             the upper triangular matrix, and the strictly lower   
00076             triangular part of A is not referenced.  If UPLO = 'L', the   
00077             leading N-by-N lower triangular part of the array A contains   
00078             the lower triangular matrix, and the strictly upper   
00079             triangular part of A is not referenced.  If DIAG = 'U', the   
00080             diagonal elements of A are also not referenced and are   
00081             assumed to be 1.   
00082             On exit, the (triangular) inverse of the original matrix, in   
00083             the same storage format.   
00084 
00085     LDA     (input) INTEGER   
00086             The leading dimension of the array A.  LDA >= max(1,N).   
00087 
00088     INFO    (output) INTEGER   
00089             = 0: successful exit   
00090             < 0: if INFO = -i, the i-th argument had an illegal value   
00091             > 0: if INFO = i, A(i,i) is exactly zero.  The triangular   
00092                  matrix is singular and its inverse can not be computed.   
00093 
00094     =====================================================================   
00095 
00096 
00097        Test the input parameters.   
00098 
00099        Parameter adjustments */
00100     /* Table of constant values */
00101      integer c__1 = 1;
00102      integer c_n1 = -1;
00103      integer c__2 = 2;
00104      Treal c_b18 = 1.;
00105      Treal c_b22 = -1.;
00106     
00107     /* System generated locals */
00108     address a__1[2];
00109     integer a_dim1, a_offset, i__1, i__2[2], i__3, i__4, i__5;
00110     char ch__1[2];
00111     /* Local variables */
00112      integer j;
00113      logical upper;
00114      integer jb, nb, nn;
00115      logical nounit;
00116 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00117 
00118 
00119     a_dim1 = *lda;
00120     a_offset = 1 + a_dim1 * 1;
00121     a -= a_offset;
00122 
00123     /* Function Body */
00124     *info = 0;
00125     upper = template_blas_lsame(uplo, "U");
00126     nounit = template_blas_lsame(diag, "N");
00127     if (! upper && ! template_blas_lsame(uplo, "L")) {
00128         *info = -1;
00129     } else if (! nounit && ! template_blas_lsame(diag, "U")) {
00130         *info = -2;
00131     } else if (*n < 0) {
00132         *info = -3;
00133     } else if (*lda < maxMACRO(1,*n)) {
00134         *info = -5;
00135     }
00136     if (*info != 0) {
00137         i__1 = -(*info);
00138         template_blas_erbla("TRTRI ", &i__1);
00139         return 0;
00140     }
00141 
00142 /*     Quick return if possible */
00143 
00144     if (*n == 0) {
00145         return 0;
00146     }
00147 
00148 /*     Check for singularity if non-unit. */
00149 
00150     if (nounit) {
00151         i__1 = *n;
00152         for (*info = 1; *info <= i__1; ++(*info)) {
00153             if (a_ref(*info, *info) == 0.) {
00154                 return 0;
00155             }
00156 /* L10: */
00157         }
00158         *info = 0;
00159     }
00160 
00161 /*     Determine the block size for this environment.   
00162 
00163    Writing concatenation */
00164     i__2[0] = 1, a__1[0] = (char*)uplo;
00165     i__2[1] = 1, a__1[1] = (char*)diag;
00166     template_blas_s_cat(ch__1, a__1, i__2, &c__2, (ftnlen)2);
00167     nb = template_lapack_ilaenv(&c__1, "DTRTRI", ch__1, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
00168             ftnlen)2);
00169     if (nb <= 1 || nb >= *n) {
00170 
00171 /*        Use unblocked code */
00172 
00173         template_lapack_trti2(uplo, diag, n, &a[a_offset], lda, info);
00174     } else {
00175 
00176 /*        Use blocked code */
00177 
00178         if (upper) {
00179 
00180 /*           Compute inverse of upper triangular matrix */
00181 
00182             i__1 = *n;
00183             i__3 = nb;
00184             for (j = 1; i__3 < 0 ? j >= i__1 : j <= i__1; j += i__3) {
00185 /* Computing MIN */
00186                 i__4 = nb, i__5 = *n - j + 1;
00187                 jb = minMACRO(i__4,i__5);
00188 
00189 /*              Compute rows 1:j-1 of current block column */
00190 
00191                 i__4 = j - 1;
00192                 template_blas_trmm("Left", "Upper", "No transpose", diag, &i__4, &jb, &
00193                         c_b18, &a[a_offset], lda, &a_ref(1, j), lda);
00194                 i__4 = j - 1;
00195                 template_blas_trsm("Right", "Upper", "No transpose", diag, &i__4, &jb, &
00196                         c_b22, &a_ref(j, j), lda, &a_ref(1, j), lda);
00197 
00198 /*              Compute inverse of current diagonal block */
00199 
00200                 template_lapack_trti2("Upper", diag, &jb, &a_ref(j, j), lda, info);
00201 /* L20: */
00202             }
00203         } else {
00204 
00205 /*           Compute inverse of lower triangular matrix */
00206 
00207             nn = (*n - 1) / nb * nb + 1;
00208             i__3 = -nb;
00209             for (j = nn; i__3 < 0 ? j >= 1 : j <= 1; j += i__3) {
00210 /* Computing MIN */
00211                 i__1 = nb, i__4 = *n - j + 1;
00212                 jb = minMACRO(i__1,i__4);
00213                 if (j + jb <= *n) {
00214 
00215 /*                 Compute rows j+jb:n of current block column */
00216 
00217                     i__1 = *n - j - jb + 1;
00218                     template_blas_trmm("Left", "Lower", "No transpose", diag, &i__1, &jb, 
00219                             &c_b18, &a_ref(j + jb, j + jb), lda, &a_ref(j + 
00220                             jb, j), lda);
00221                     i__1 = *n - j - jb + 1;
00222                     template_blas_trsm("Right", "Lower", "No transpose", diag, &i__1, &jb,
00223                              &c_b22, &a_ref(j, j), lda, &a_ref(j + jb, j), 
00224                             lda);
00225                 }
00226 
00227 /*              Compute inverse of current diagonal block */
00228 
00229                 template_lapack_trti2("Lower", diag, &jb, &a_ref(j, j), lda, info);
00230 /* L30: */
00231             }
00232         }
00233     }
00234 
00235     return 0;
00236 
00237 /*     End of DTRTRI */
00238 
00239 } /* dtrtri_ */
00240 
00241 #undef a_ref
00242 
00243 
00244 #endif

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1