ergo
|
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