ergo
template_lapack_larra.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_LARRA_HEADER
00036 #define TEMPLATE_LAPACK_LARRA_HEADER
00037 
00038 template<class Treal>
00039 int template_lapack_larra(const integer *n, Treal *d__, Treal *e, 
00040         Treal *e2, Treal *spltol, Treal *tnrm, integer *nsplit, 
00041          integer *isplit, integer *info)
00042 {
00043     /* System generated locals */
00044     integer i__1;
00045     Treal d__1, d__2;
00046 
00047 
00048     /* Local variables */
00049     integer i__;
00050     Treal tmp1, eabs;
00051 
00052 
00053 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00054 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00055 /*     November 2006 */
00056 
00057 /*     .. Scalar Arguments .. */
00058 /*     .. */
00059 /*     .. Array Arguments .. */
00060 /*     .. */
00061 
00062 /*  Purpose */
00063 /*  ======= */
00064 
00065 /*  Compute the splitting points with threshold SPLTOL. */
00066 /*  DLARRA sets any "small" off-diagonal elements to zero. */
00067 
00068 /*  Arguments */
00069 /*  ========= */
00070 
00071 /*  N       (input) INTEGER */
00072 /*          The order of the matrix. N > 0. */
00073 
00074 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
00075 /*          On entry, the N diagonal elements of the tridiagonal */
00076 /*          matrix T. */
00077 
00078 /*  E       (input/output) DOUBLE PRECISION array, dimension (N) */
00079 /*          On entry, the first (N-1) entries contain the subdiagonal */
00080 /*          elements of the tridiagonal matrix T; E(N) need not be set. */
00081 /*          On exit, the entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT, */
00082 /*          are set to zero, the other entries of E are untouched. */
00083 
00084 /*  E2      (input/output) DOUBLE PRECISION array, dimension (N) */
00085 /*          On entry, the first (N-1) entries contain the SQUARES of the */
00086 /*          subdiagonal elements of the tridiagonal matrix T; */
00087 /*          E2(N) need not be set. */
00088 /*          On exit, the entries E2( ISPLIT( I ) ), */
00089 /*          1 <= I <= NSPLIT, have been set to zero */
00090 
00091 /*  SPLTOL (input) DOUBLE PRECISION */
00092 /*          The threshold for splitting. Two criteria can be used: */
00093 /*          SPLTOL<0 : criterion based on absolute off-diagonal value */
00094 /*          SPLTOL>0 : criterion that preserves relative accuracy */
00095 
00096 /*  TNRM (input) DOUBLE PRECISION */
00097 /*          The norm of the matrix. */
00098 
00099 /*  NSPLIT  (output) INTEGER */
00100 /*          The number of blocks T splits into. 1 <= NSPLIT <= N. */
00101 
00102 /*  ISPLIT  (output) INTEGER array, dimension (N) */
00103 /*          The splitting points, at which T breaks up into blocks. */
00104 /*          The first block consists of rows/columns 1 to ISPLIT(1), */
00105 /*          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
00106 /*          etc., and the NSPLIT-th consists of rows/columns */
00107 /*          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
00108 
00109 
00110 /*  INFO    (output) INTEGER */
00111 /*          = 0:  successful exit */
00112 
00113 /*  Further Details */
00114 /*  =============== */
00115 
00116 /*  Based on contributions by */
00117 /*     Beresford Parlett, University of California, Berkeley, USA */
00118 /*     Jim Demmel, University of California, Berkeley, USA */
00119 /*     Inderjit Dhillon, University of Texas, Austin, USA */
00120 /*     Osni Marques, LBNL/NERSC, USA */
00121 /*     Christof Voemel, University of California, Berkeley, USA */
00122 
00123 /*  ===================================================================== */
00124 
00125 /*     .. Parameters .. */
00126 /*     .. */
00127 /*     .. Local Scalars .. */
00128 /*     .. */
00129 /*     .. Intrinsic Functions .. */
00130 /*     .. */
00131 /*     .. Executable Statements .. */
00132 
00133     /* Parameter adjustments */
00134     --isplit;
00135     --e2;
00136     --e;
00137     --d__;
00138 
00139     /* Function Body */
00140     *info = 0;
00141 /*     Compute splitting points */
00142     *nsplit = 1;
00143     if (*spltol < 0.) {
00144 /*        Criterion based on absolute off-diagonal value */
00145         tmp1 = absMACRO(*spltol) * *tnrm;
00146         i__1 = *n - 1;
00147         for (i__ = 1; i__ <= i__1; ++i__) {
00148             eabs = (d__1 = e[i__], absMACRO(d__1));
00149             if (eabs <= tmp1) {
00150                 e[i__] = 0.;
00151                 e2[i__] = 0.;
00152                 isplit[*nsplit] = i__;
00153                 ++(*nsplit);
00154             }
00155 /* L9: */
00156         }
00157     } else {
00158 /*        Criterion that guarantees relative accuracy */
00159         i__1 = *n - 1;
00160         for (i__ = 1; i__ <= i__1; ++i__) {
00161             eabs = (d__1 = e[i__], absMACRO(d__1));
00162             if (eabs <= *spltol * template_blas_sqrt((d__1 = d__[i__], absMACRO(d__1))) * template_blas_sqrt((
00163                     d__2 = d__[i__ + 1], absMACRO(d__2)))) {
00164                 e[i__] = 0.;
00165                 e2[i__] = 0.;
00166                 isplit[*nsplit] = i__;
00167                 ++(*nsplit);
00168             }
00169 /* L10: */
00170         }
00171     }
00172     isplit[*nsplit] = *n;
00173     return 0;
00174 
00175 /*     End of DLARRA */
00176 
00177 } /* dlarra_ */
00178 
00179 #endif