ergo
template_lapack_lagts.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_LAGTS_HEADER
00036 #define TEMPLATE_LAPACK_LAGTS_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_lagts(const integer *job, const integer *n, const Treal *a, 
00041         const Treal *b, const Treal *c__, const Treal *d__, const integer *in, 
00042         Treal *y, Treal *tol, integer *info)
00043 {
00044 /*  -- LAPACK auxiliary routine (version 3.0) --   
00045        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00046        Courant Institute, Argonne National Lab, and Rice University   
00047        October 31, 1992   
00048 
00049 
00050     Purpose   
00051     =======   
00052 
00053     DLAGTS may be used to solve one of the systems of equations   
00054 
00055        (T - lambda*I)*x = y   or   (T - lambda*I)'*x = y,   
00056 
00057     where T is an n by n tridiagonal matrix, for x, following the   
00058     factorization of (T - lambda*I) as   
00059 
00060        (T - lambda*I) = P*L*U ,   
00061 
00062     by routine DLAGTF. The choice of equation to be solved is   
00063     controlled by the argument JOB, and in each case there is an option   
00064     to perturb zero or very small diagonal elements of U, this option   
00065     being intended for use in applications such as inverse iteration.   
00066 
00067     Arguments   
00068     =========   
00069 
00070     JOB     (input) INTEGER   
00071             Specifies the job to be performed by DLAGTS as follows:   
00072             =  1: The equations  (T - lambda*I)x = y  are to be solved,   
00073                   but diagonal elements of U are not to be perturbed.   
00074             = -1: The equations  (T - lambda*I)x = y  are to be solved   
00075                   and, if overflow would otherwise occur, the diagonal   
00076                   elements of U are to be perturbed. See argument TOL   
00077                   below.   
00078             =  2: The equations  (T - lambda*I)'x = y  are to be solved,   
00079                   but diagonal elements of U are not to be perturbed.   
00080             = -2: The equations  (T - lambda*I)'x = y  are to be solved   
00081                   and, if overflow would otherwise occur, the diagonal   
00082                   elements of U are to be perturbed. See argument TOL   
00083                   below.   
00084 
00085     N       (input) INTEGER   
00086             The order of the matrix T.   
00087 
00088     A       (input) DOUBLE PRECISION array, dimension (N)   
00089             On entry, A must contain the diagonal elements of U as   
00090             returned from DLAGTF.   
00091 
00092     B       (input) DOUBLE PRECISION array, dimension (N-1)   
00093             On entry, B must contain the first super-diagonal elements of   
00094             U as returned from DLAGTF.   
00095 
00096     C       (input) DOUBLE PRECISION array, dimension (N-1)   
00097             On entry, C must contain the sub-diagonal elements of L as   
00098             returned from DLAGTF.   
00099 
00100     D       (input) DOUBLE PRECISION array, dimension (N-2)   
00101             On entry, D must contain the second super-diagonal elements   
00102             of U as returned from DLAGTF.   
00103 
00104     IN      (input) INTEGER array, dimension (N)   
00105             On entry, IN must contain details of the matrix P as returned   
00106             from DLAGTF.   
00107 
00108     Y       (input/output) DOUBLE PRECISION array, dimension (N)   
00109             On entry, the right hand side vector y.   
00110             On exit, Y is overwritten by the solution vector x.   
00111 
00112     TOL     (input/output) DOUBLE PRECISION   
00113             On entry, with  JOB .lt. 0, TOL should be the minimum   
00114             perturbation to be made to very small diagonal elements of U.   
00115             TOL should normally be chosen as about eps*norm(U), where eps   
00116             is the relative machine precision, but if TOL is supplied as   
00117             non-positive, then it is reset to eps*max( abs( u(i,j) ) ).   
00118             If  JOB .gt. 0  then TOL is not referenced.   
00119 
00120             On exit, TOL is changed as described above, only if TOL is   
00121             non-positive on entry. Otherwise TOL is unchanged.   
00122 
00123     INFO    (output) INTEGER   
00124             = 0   : successful exit   
00125             .lt. 0: if INFO = -i, the i-th argument had an illegal value   
00126             .gt. 0: overflow would occur when computing the INFO(th)   
00127                     element of the solution vector x. This can only occur   
00128                     when JOB is supplied as positive and either means   
00129                     that a diagonal element of U is very small, or that   
00130                     the elements of the right-hand side vector y are very   
00131                     large.   
00132 
00133     =====================================================================   
00134 
00135 
00136        Parameter adjustments */
00137     /* System generated locals */
00138     integer i__1;
00139     Treal d__1, d__2, d__3, d__4, d__5;
00140     /* Local variables */
00141      Treal temp, pert;
00142      integer k;
00143      Treal absak, sfmin, ak;
00144      Treal bignum, eps;
00145 
00146     --y;
00147     --in;
00148     --d__;
00149     --c__;
00150     --b;
00151     --a;
00152 
00153     /* Function Body */
00154     *info = 0;
00155     if (absMACRO(*job) > 2 || *job == 0) {
00156         *info = -1;
00157     } else if (*n < 0) {
00158         *info = -2;
00159     }
00160     if (*info != 0) {
00161         i__1 = -(*info);
00162         template_blas_erbla("LAGTS ", &i__1);
00163         return 0;
00164     }
00165 
00166     if (*n == 0) {
00167         return 0;
00168     }
00169 
00170     eps = template_lapack_lamch("Epsilon", (Treal)0);
00171     sfmin = template_lapack_lamch("Safe minimum", (Treal)0);
00172     bignum = 1. / sfmin;
00173 
00174     if (*job < 0) {
00175         if (*tol <= 0.) {
00176             *tol = absMACRO(a[1]);
00177             if (*n > 1) {
00178 /* Computing MAX */
00179                 d__1 = *tol, d__2 = absMACRO(a[2]), d__1 = maxMACRO(d__1,d__2), d__2 = 
00180                         absMACRO(b[1]);
00181                 *tol = maxMACRO(d__1,d__2);
00182             }
00183             i__1 = *n;
00184             for (k = 3; k <= i__1; ++k) {
00185 /* Computing MAX */
00186                 d__4 = *tol, d__5 = (d__1 = a[k], absMACRO(d__1)), d__4 = maxMACRO(d__4,
00187                         d__5), d__5 = (d__2 = b[k - 1], absMACRO(d__2)), d__4 = 
00188                         maxMACRO(d__4,d__5), d__5 = (d__3 = d__[k - 2], absMACRO(d__3));
00189                 *tol = maxMACRO(d__4,d__5);
00190 /* L10: */
00191             }
00192             *tol *= eps;
00193             if (*tol == 0.) {
00194                 *tol = eps;
00195             }
00196         }
00197     }
00198 
00199     if (absMACRO(*job) == 1) {
00200         i__1 = *n;
00201         for (k = 2; k <= i__1; ++k) {
00202             if (in[k - 1] == 0) {
00203                 y[k] -= c__[k - 1] * y[k - 1];
00204             } else {
00205                 temp = y[k - 1];
00206                 y[k - 1] = y[k];
00207                 y[k] = temp - c__[k - 1] * y[k];
00208             }
00209 /* L20: */
00210         }
00211         if (*job == 1) {
00212             for (k = *n; k >= 1; --k) {
00213                 if (k <= *n - 2) {
00214                     temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
00215                 } else if (k == *n - 1) {
00216                     temp = y[k] - b[k] * y[k + 1];
00217                 } else {
00218                     temp = y[k];
00219                 }
00220                 ak = a[k];
00221                 absak = absMACRO(ak);
00222                 if (absak < 1.) {
00223                     if (absak < sfmin) {
00224                         if (absak == 0. || absMACRO(temp) * sfmin > absak) {
00225                             *info = k;
00226                             return 0;
00227                         } else {
00228                             temp *= bignum;
00229                             ak *= bignum;
00230                         }
00231                     } else if (absMACRO(temp) > absak * bignum) {
00232                         *info = k;
00233                         return 0;
00234                     }
00235                 }
00236                 y[k] = temp / ak;
00237 /* L30: */
00238             }
00239         } else {
00240             for (k = *n; k >= 1; --k) {
00241                 if (k <= *n - 2) {
00242                     temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
00243                 } else if (k == *n - 1) {
00244                     temp = y[k] - b[k] * y[k + 1];
00245                 } else {
00246                     temp = y[k];
00247                 }
00248                 ak = a[k];
00249                 pert = template_lapack_d_sign(tol, &ak);
00250 L40:
00251                 absak = absMACRO(ak);
00252                 if (absak < 1.) {
00253                     if (absak < sfmin) {
00254                         if (absak == 0. || absMACRO(temp) * sfmin > absak) {
00255                             ak += pert;
00256                             pert *= 2;
00257                             goto L40;
00258                         } else {
00259                             temp *= bignum;
00260                             ak *= bignum;
00261                         }
00262                     } else if (absMACRO(temp) > absak * bignum) {
00263                         ak += pert;
00264                         pert *= 2;
00265                         goto L40;
00266                     }
00267                 }
00268                 y[k] = temp / ak;
00269 /* L50: */
00270             }
00271         }
00272     } else {
00273 
00274 /*        Come to here if  JOB = 2 or -2 */
00275 
00276         if (*job == 2) {
00277             i__1 = *n;
00278             for (k = 1; k <= i__1; ++k) {
00279                 if (k >= 3) {
00280                     temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
00281                 } else if (k == 2) {
00282                     temp = y[k] - b[k - 1] * y[k - 1];
00283                 } else {
00284                     temp = y[k];
00285                 }
00286                 ak = a[k];
00287                 absak = absMACRO(ak);
00288                 if (absak < 1.) {
00289                     if (absak < sfmin) {
00290                         if (absak == 0. || absMACRO(temp) * sfmin > absak) {
00291                             *info = k;
00292                             return 0;
00293                         } else {
00294                             temp *= bignum;
00295                             ak *= bignum;
00296                         }
00297                     } else if (absMACRO(temp) > absak * bignum) {
00298                         *info = k;
00299                         return 0;
00300                     }
00301                 }
00302                 y[k] = temp / ak;
00303 /* L60: */
00304             }
00305         } else {
00306             i__1 = *n;
00307             for (k = 1; k <= i__1; ++k) {
00308                 if (k >= 3) {
00309                     temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
00310                 } else if (k == 2) {
00311                     temp = y[k] - b[k - 1] * y[k - 1];
00312                 } else {
00313                     temp = y[k];
00314                 }
00315                 ak = a[k];
00316                 pert = template_lapack_d_sign(tol, &ak);
00317 L70:
00318                 absak = absMACRO(ak);
00319                 if (absak < 1.) {
00320                     if (absak < sfmin) {
00321                         if (absak == 0. || absMACRO(temp) * sfmin > absak) {
00322                             ak += pert;
00323                             pert *= 2;
00324                             goto L70;
00325                         } else {
00326                             temp *= bignum;
00327                             ak *= bignum;
00328                         }
00329                     } else if (absMACRO(temp) > absak * bignum) {
00330                         ak += pert;
00331                         pert *= 2;
00332                         goto L70;
00333                     }
00334                 }
00335                 y[k] = temp / ak;
00336 /* L80: */
00337             }
00338         }
00339 
00340         for (k = *n; k >= 2; --k) {
00341             if (in[k - 1] == 0) {
00342                 y[k - 1] -= c__[k - 1] * y[k];
00343             } else {
00344                 temp = y[k - 1];
00345                 y[k - 1] = y[k];
00346                 y[k] = temp - c__[k - 1] * y[k];
00347             }
00348 /* L90: */
00349         }
00350     }
00351 
00352 /*     End of DLAGTS */
00353 
00354     return 0;
00355 } /* dlagts_ */
00356 
00357 #endif