template_lapack_lagtf.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_LAGTF_HEADER
00036 #define TEMPLATE_LAPACK_LAGTF_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_lagtf(const integer *n, Treal *a, const Treal *lambda, 
00041         Treal *b, Treal *c__, const Treal *tol, Treal *d__, 
00042         integer *in, 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        June 30, 1999   
00048 
00049 
00050     Purpose   
00051     =======   
00052 
00053     DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n   
00054     tridiagonal matrix and lambda is a scalar, as   
00055 
00056        T - lambda*I = PLU,   
00057 
00058     where P is a permutation matrix, L is a unit lower tridiagonal matrix   
00059     with at most one non-zero sub-diagonal elements per column and U is   
00060     an upper triangular matrix with at most two non-zero super-diagonal   
00061     elements per column.   
00062 
00063     The factorization is obtained by Gaussian elimination with partial   
00064     pivoting and implicit row scaling.   
00065 
00066     The parameter LAMBDA is included in the routine so that DLAGTF may   
00067     be used, in conjunction with DLAGTS, to obtain eigenvectors of T by   
00068     inverse iteration.   
00069 
00070     Arguments   
00071     =========   
00072 
00073     N       (input) INTEGER   
00074             The order of the matrix T.   
00075 
00076     A       (input/output) DOUBLE PRECISION array, dimension (N)   
00077             On entry, A must contain the diagonal elements of T.   
00078 
00079             On exit, A is overwritten by the n diagonal elements of the   
00080             upper triangular matrix U of the factorization of T.   
00081 
00082     LAMBDA  (input) DOUBLE PRECISION   
00083             On entry, the scalar lambda.   
00084 
00085     B       (input/output) DOUBLE PRECISION array, dimension (N-1)   
00086             On entry, B must contain the (n-1) super-diagonal elements of   
00087             T.   
00088 
00089             On exit, B is overwritten by the (n-1) super-diagonal   
00090             elements of the matrix U of the factorization of T.   
00091 
00092     C       (input/output) DOUBLE PRECISION array, dimension (N-1)   
00093             On entry, C must contain the (n-1) sub-diagonal elements of   
00094             T.   
00095 
00096             On exit, C is overwritten by the (n-1) sub-diagonal elements   
00097             of the matrix L of the factorization of T.   
00098 
00099     TOL     (input) DOUBLE PRECISION   
00100             On entry, a relative tolerance used to indicate whether or   
00101             not the matrix (T - lambda*I) is nearly singular. TOL should   
00102             normally be chose as approximately the largest relative error   
00103             in the elements of T. For example, if the elements of T are   
00104             correct to about 4 significant figures, then TOL should be   
00105             set to about 5*10**(-4). If TOL is supplied as less than eps,   
00106             where eps is the relative machine precision, then the value   
00107             eps is used in place of TOL.   
00108 
00109     D       (output) DOUBLE PRECISION array, dimension (N-2)   
00110             On exit, D is overwritten by the (n-2) second super-diagonal   
00111             elements of the matrix U of the factorization of T.   
00112 
00113     IN      (output) INTEGER array, dimension (N)   
00114             On exit, IN contains details of the permutation matrix P. If   
00115             an interchange occurred at the kth step of the elimination,   
00116             then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)   
00117             returns the smallest positive integer j such that   
00118 
00119                abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,   
00120 
00121             where norm( A(j) ) denotes the sum of the absolute values of   
00122             the jth row of the matrix A. If no such j exists then IN(n)   
00123             is returned as zero. If IN(n) is returned as positive, then a   
00124             diagonal element of U is small, indicating that   
00125             (T - lambda*I) is singular or nearly singular,   
00126 
00127     INFO    (output) INTEGER   
00128             = 0   : successful exit   
00129             .lt. 0: if INFO = -k, the kth argument had an illegal value   
00130 
00131    =====================================================================   
00132 
00133 
00134        Parameter adjustments */
00135     /* System generated locals */
00136     integer i__1;
00137     Treal d__1, d__2;
00138     /* Local variables */
00139      Treal temp, mult;
00140      integer k;
00141      Treal scale1, scale2;
00142      Treal tl;
00143      Treal eps, piv1, piv2;
00144 
00145     --in;
00146     --d__;
00147     --c__;
00148     --b;
00149     --a;
00150 
00151     /* Function Body */
00152     *info = 0;
00153     if (*n < 0) {
00154         *info = -1;
00155         i__1 = -(*info);
00156         template_blas_erbla("LAGTF ", &i__1);
00157         return 0;
00158     }
00159 
00160     if (*n == 0) {
00161         return 0;
00162     }
00163 
00164     a[1] -= *lambda;
00165     in[*n] = 0;
00166     if (*n == 1) {
00167         if (a[1] == 0.) {
00168             in[1] = 1;
00169         }
00170         return 0;
00171     }
00172 
00173     eps = template_lapack_lamch("Epsilon", (Treal)0);
00174 
00175     tl = maxMACRO(*tol,eps);
00176     scale1 = absMACRO(a[1]) + absMACRO(b[1]);
00177     i__1 = *n - 1;
00178     for (k = 1; k <= i__1; ++k) {
00179         a[k + 1] -= *lambda;
00180         scale2 = (d__1 = c__[k], absMACRO(d__1)) + (d__2 = a[k + 1], absMACRO(d__2));
00181         if (k < *n - 1) {
00182             scale2 += (d__1 = b[k + 1], absMACRO(d__1));
00183         }
00184         if (a[k] == 0.) {
00185             piv1 = 0.;
00186         } else {
00187             piv1 = (d__1 = a[k], absMACRO(d__1)) / scale1;
00188         }
00189         if (c__[k] == 0.) {
00190             in[k] = 0;
00191             piv2 = 0.;
00192             scale1 = scale2;
00193             if (k < *n - 1) {
00194                 d__[k] = 0.;
00195             }
00196         } else {
00197             piv2 = (d__1 = c__[k], absMACRO(d__1)) / scale2;
00198             if (piv2 <= piv1) {
00199                 in[k] = 0;
00200                 scale1 = scale2;
00201                 c__[k] /= a[k];
00202                 a[k + 1] -= c__[k] * b[k];
00203                 if (k < *n - 1) {
00204                     d__[k] = 0.;
00205                 }
00206             } else {
00207                 in[k] = 1;
00208                 mult = a[k] / c__[k];
00209                 a[k] = c__[k];
00210                 temp = a[k + 1];
00211                 a[k + 1] = b[k] - mult * temp;
00212                 if (k < *n - 1) {
00213                     d__[k] = b[k + 1];
00214                     b[k + 1] = -mult * d__[k];
00215                 }
00216                 b[k] = temp;
00217                 c__[k] = mult;
00218             }
00219         }
00220         if (maxMACRO(piv1,piv2) <= tl && in[*n] == 0) {
00221             in[*n] = k;
00222         }
00223 /* L10: */
00224     }
00225     if ((d__1 = a[*n], absMACRO(d__1)) <= scale1 * tl && in[*n] == 0) {
00226         in[*n] = *n;
00227     }
00228 
00229     return 0;
00230 
00231 /*     End of DLAGTF */
00232 
00233 } /* dlagtf_ */
00234 
00235 #endif

Generated on Mon Sep 17 14:32:56 2012 for ergo by  doxygen 1.4.7