00001 /* Ergo, version 3.7, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 00004 * and Anastasia Kruchinina. 00005 * 00006 * This program is free software: you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation, either version 3 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00018 * 00019 * Primary academic reference: 00020 * Ergo: An open-source program for linear-scaling electronic structure 00021 * calculations, 00022 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 00023 * Kruchinina, 00024 * SoftwareX 7, 107 (2018), 00025 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 00026 * 00027 * For further information about Ergo, see <http://www.ergoscf.org>. 00028 */ 00029 00030 /* This file belongs to the template_lapack part of the Ergo source 00031 * code. The source files in the template_lapack directory are modified 00032 * versions of files originally distributed as CLAPACK, see the 00033 * Copyright/license notice in the file template_lapack/COPYING. 00034 */ 00035 00036 00037 #ifndef TEMPLATE_LAPACK_LARRA_HEADER 00038 #define TEMPLATE_LAPACK_LARRA_HEADER 00039 00040 template<class Treal> 00041 int template_lapack_larra(const integer *n, Treal *d__, Treal *e, 00042 Treal *e2, Treal *spltol, Treal *tnrm, integer *nsplit, 00043 integer *isplit, integer *info) 00044 { 00045 /* System generated locals */ 00046 integer i__1; 00047 Treal d__1, d__2; 00048 00049 00050 /* Local variables */ 00051 integer i__; 00052 Treal tmp1, eabs; 00053 00054 00055 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00056 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00057 /* November 2006 */ 00058 00059 /* .. Scalar Arguments .. */ 00060 /* .. */ 00061 /* .. Array Arguments .. */ 00062 /* .. */ 00063 00064 /* Purpose */ 00065 /* ======= */ 00066 00067 /* Compute the splitting points with threshold SPLTOL. */ 00068 /* DLARRA sets any "small" off-diagonal elements to zero. */ 00069 00070 /* Arguments */ 00071 /* ========= */ 00072 00073 /* N (input) INTEGER */ 00074 /* The order of the matrix. N > 0. */ 00075 00076 /* D (input) DOUBLE PRECISION array, dimension (N) */ 00077 /* On entry, the N diagonal elements of the tridiagonal */ 00078 /* matrix T. */ 00079 00080 /* E (input/output) DOUBLE PRECISION array, dimension (N) */ 00081 /* On entry, the first (N-1) entries contain the subdiagonal */ 00082 /* elements of the tridiagonal matrix T; E(N) need not be set. */ 00083 /* On exit, the entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT, */ 00084 /* are set to zero, the other entries of E are untouched. */ 00085 00086 /* E2 (input/output) DOUBLE PRECISION array, dimension (N) */ 00087 /* On entry, the first (N-1) entries contain the SQUARES of the */ 00088 /* subdiagonal elements of the tridiagonal matrix T; */ 00089 /* E2(N) need not be set. */ 00090 /* On exit, the entries E2( ISPLIT( I ) ), */ 00091 /* 1 <= I <= NSPLIT, have been set to zero */ 00092 00093 /* SPLTOL (input) DOUBLE PRECISION */ 00094 /* The threshold for splitting. Two criteria can be used: */ 00095 /* SPLTOL<0 : criterion based on absolute off-diagonal value */ 00096 /* SPLTOL>0 : criterion that preserves relative accuracy */ 00097 00098 /* TNRM (input) DOUBLE PRECISION */ 00099 /* The norm of the matrix. */ 00100 00101 /* NSPLIT (output) INTEGER */ 00102 /* The number of blocks T splits into. 1 <= NSPLIT <= N. */ 00103 00104 /* ISPLIT (output) INTEGER array, dimension (N) */ 00105 /* The splitting points, at which T breaks up into blocks. */ 00106 /* The first block consists of rows/columns 1 to ISPLIT(1), */ 00107 /* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */ 00108 /* etc., and the NSPLIT-th consists of rows/columns */ 00109 /* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */ 00110 00111 00112 /* INFO (output) INTEGER */ 00113 /* = 0: successful exit */ 00114 00115 /* Further Details */ 00116 /* =============== */ 00117 00118 /* Based on contributions by */ 00119 /* Beresford Parlett, University of California, Berkeley, USA */ 00120 /* Jim Demmel, University of California, Berkeley, USA */ 00121 /* Inderjit Dhillon, University of Texas, Austin, USA */ 00122 /* Osni Marques, LBNL/NERSC, USA */ 00123 /* Christof Voemel, University of California, Berkeley, USA */ 00124 00125 /* ===================================================================== */ 00126 00127 /* .. Parameters .. */ 00128 /* .. */ 00129 /* .. Local Scalars .. */ 00130 /* .. */ 00131 /* .. Intrinsic Functions .. */ 00132 /* .. */ 00133 /* .. Executable Statements .. */ 00134 00135 /* Parameter adjustments */ 00136 --isplit; 00137 --e2; 00138 --e; 00139 --d__; 00140 00141 /* Function Body */ 00142 *info = 0; 00143 /* Compute splitting points */ 00144 *nsplit = 1; 00145 if (*spltol < 0.) { 00146 /* Criterion based on absolute off-diagonal value */ 00147 tmp1 = absMACRO(*spltol) * *tnrm; 00148 i__1 = *n - 1; 00149 for (i__ = 1; i__ <= i__1; ++i__) { 00150 eabs = (d__1 = e[i__], absMACRO(d__1)); 00151 if (eabs <= tmp1) { 00152 e[i__] = 0.; 00153 e2[i__] = 0.; 00154 isplit[*nsplit] = i__; 00155 ++(*nsplit); 00156 } 00157 /* L9: */ 00158 } 00159 } else { 00160 /* Criterion that guarantees relative accuracy */ 00161 i__1 = *n - 1; 00162 for (i__ = 1; i__ <= i__1; ++i__) { 00163 eabs = (d__1 = e[i__], absMACRO(d__1)); 00164 if (eabs <= *spltol * template_blas_sqrt((d__1 = d__[i__], absMACRO(d__1))) * template_blas_sqrt(( 00165 d__2 = d__[i__ + 1], absMACRO(d__2)))) { 00166 e[i__] = 0.; 00167 e2[i__] = 0.; 00168 isplit[*nsplit] = i__; 00169 ++(*nsplit); 00170 } 00171 /* L10: */ 00172 } 00173 } 00174 isplit[*nsplit] = *n; 00175 return 0; 00176 00177 /* End of DLARRA */ 00178 00179 } /* dlarra_ */ 00180 00181 #endif