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_LASCL_HEADER 00036 #define TEMPLATE_LAPACK_LASCL_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_lascl(const char *type__, const integer *kl, const integer *ku, 00041 const Treal *cfrom, const Treal *cto, const integer *m, const integer *n, 00042 Treal *a, const integer *lda, 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 February 29, 1992 00048 00049 00050 Purpose 00051 ======= 00052 00053 DLASCL multiplies the M by N real matrix A by the real scalar 00054 CTO/CFROM. This is done without over/underflow as long as the final 00055 result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that 00056 A may be full, upper triangular, lower triangular, upper Hessenberg, 00057 or banded. 00058 00059 Arguments 00060 ========= 00061 00062 TYPE (input) CHARACTER*1 00063 TYPE indices the storage type of the input matrix. 00064 = 'G': A is a full matrix. 00065 = 'L': A is a lower triangular matrix. 00066 = 'U': A is an upper triangular matrix. 00067 = 'H': A is an upper Hessenberg matrix. 00068 = 'B': A is a symmetric band matrix with lower bandwidth KL 00069 and upper bandwidth KU and with the only the lower 00070 half stored. 00071 = 'Q': A is a symmetric band matrix with lower bandwidth KL 00072 and upper bandwidth KU and with the only the upper 00073 half stored. 00074 = 'Z': A is a band matrix with lower bandwidth KL and upper 00075 bandwidth KU. 00076 00077 KL (input) INTEGER 00078 The lower bandwidth of A. Referenced only if TYPE = 'B', 00079 'Q' or 'Z'. 00080 00081 KU (input) INTEGER 00082 The upper bandwidth of A. Referenced only if TYPE = 'B', 00083 'Q' or 'Z'. 00084 00085 CFROM (input) DOUBLE PRECISION 00086 CTO (input) DOUBLE PRECISION 00087 The matrix A is multiplied by CTO/CFROM. A(I,J) is computed 00088 without over/underflow if the final result CTO*A(I,J)/CFROM 00089 can be represented without over/underflow. CFROM must be 00090 nonzero. 00091 00092 M (input) INTEGER 00093 The number of rows of the matrix A. M >= 0. 00094 00095 N (input) INTEGER 00096 The number of columns of the matrix A. N >= 0. 00097 00098 A (input/output) DOUBLE PRECISION array, dimension (LDA,M) 00099 The matrix to be multiplied by CTO/CFROM. See TYPE for the 00100 storage type. 00101 00102 LDA (input) INTEGER 00103 The leading dimension of the array A. LDA >= max(1,M). 00104 00105 INFO (output) INTEGER 00106 0 - successful exit 00107 <0 - if INFO = -i, the i-th argument had an illegal value. 00108 00109 ===================================================================== 00110 00111 00112 Test the input arguments 00113 00114 Parameter adjustments */ 00115 /* System generated locals */ 00116 integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5; 00117 /* Local variables */ 00118 logical done; 00119 Treal ctoc; 00120 integer i__, j; 00121 integer itype, k1, k2, k3, k4; 00122 Treal cfrom1; 00123 Treal cfromc; 00124 Treal bignum, smlnum, mul, cto1; 00125 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00126 00127 a_dim1 = *lda; 00128 a_offset = 1 + a_dim1 * 1; 00129 a -= a_offset; 00130 00131 /* Function Body */ 00132 *info = 0; 00133 00134 if (template_blas_lsame(type__, "G")) { 00135 itype = 0; 00136 } else if (template_blas_lsame(type__, "L")) { 00137 itype = 1; 00138 } else if (template_blas_lsame(type__, "U")) { 00139 itype = 2; 00140 } else if (template_blas_lsame(type__, "H")) { 00141 itype = 3; 00142 } else if (template_blas_lsame(type__, "B")) { 00143 itype = 4; 00144 } else if (template_blas_lsame(type__, "Q")) { 00145 itype = 5; 00146 } else if (template_blas_lsame(type__, "Z")) { 00147 itype = 6; 00148 } else { 00149 itype = -1; 00150 } 00151 00152 if (itype == -1) { 00153 *info = -1; 00154 } else if (*cfrom == 0.) { 00155 *info = -4; 00156 } else if (*m < 0) { 00157 *info = -6; 00158 } else if (*n < 0 || ( itype == 4 && *n != *m ) || ( itype == 5 && *n != *m ) ) { 00159 *info = -7; 00160 } else if (itype <= 3 && *lda < maxMACRO(1,*m)) { 00161 *info = -9; 00162 } else if (itype >= 4) { 00163 /* Computing MAX */ 00164 i__1 = *m - 1; 00165 if (*kl < 0 || *kl > maxMACRO(i__1,0)) { 00166 *info = -2; 00167 } else /* if(complicated condition) */ { 00168 /* Computing MAX */ 00169 i__1 = *n - 1; 00170 if (*ku < 0 || *ku > maxMACRO(i__1,0) || ( (itype == 4 || itype == 5) && 00171 *kl != *ku ) ) { 00172 *info = -3; 00173 } else if ( ( itype == 4 && *lda < *kl + 1 ) || ( itype == 5 && *lda < * 00174 ku + 1 ) || ( itype == 6 && *lda < (*kl << 1) + *ku + 1 ) ) { 00175 *info = -9; 00176 } 00177 } 00178 } 00179 00180 if (*info != 0) { 00181 i__1 = -(*info); 00182 template_blas_erbla("LASCL ", &i__1); 00183 return 0; 00184 } 00185 00186 /* Quick return if possible */ 00187 00188 if (*n == 0 || *m == 0) { 00189 return 0; 00190 } 00191 00192 /* Get machine parameters */ 00193 00194 smlnum = template_lapack_lamch("S", (Treal)0); 00195 bignum = 1. / smlnum; 00196 00197 cfromc = *cfrom; 00198 ctoc = *cto; 00199 00200 L10: 00201 cfrom1 = cfromc * smlnum; 00202 cto1 = ctoc / bignum; 00203 if (absMACRO(cfrom1) > absMACRO(ctoc) && ctoc != 0.) { 00204 mul = smlnum; 00205 done = FALSE_; 00206 cfromc = cfrom1; 00207 } else if (absMACRO(cto1) > absMACRO(cfromc)) { 00208 mul = bignum; 00209 done = FALSE_; 00210 ctoc = cto1; 00211 } else { 00212 mul = ctoc / cfromc; 00213 done = TRUE_; 00214 } 00215 00216 if (itype == 0) { 00217 00218 /* Full matrix */ 00219 00220 i__1 = *n; 00221 for (j = 1; j <= i__1; ++j) { 00222 i__2 = *m; 00223 for (i__ = 1; i__ <= i__2; ++i__) { 00224 a_ref(i__, j) = a_ref(i__, j) * mul; 00225 /* L20: */ 00226 } 00227 /* L30: */ 00228 } 00229 00230 } else if (itype == 1) { 00231 00232 /* Lower triangular matrix */ 00233 00234 i__1 = *n; 00235 for (j = 1; j <= i__1; ++j) { 00236 i__2 = *m; 00237 for (i__ = j; i__ <= i__2; ++i__) { 00238 a_ref(i__, j) = a_ref(i__, j) * mul; 00239 /* L40: */ 00240 } 00241 /* L50: */ 00242 } 00243 00244 } else if (itype == 2) { 00245 00246 /* Upper triangular matrix */ 00247 00248 i__1 = *n; 00249 for (j = 1; j <= i__1; ++j) { 00250 i__2 = minMACRO(j,*m); 00251 for (i__ = 1; i__ <= i__2; ++i__) { 00252 a_ref(i__, j) = a_ref(i__, j) * mul; 00253 /* L60: */ 00254 } 00255 /* L70: */ 00256 } 00257 00258 } else if (itype == 3) { 00259 00260 /* Upper Hessenberg matrix */ 00261 00262 i__1 = *n; 00263 for (j = 1; j <= i__1; ++j) { 00264 /* Computing MIN */ 00265 i__3 = j + 1; 00266 i__2 = minMACRO(i__3,*m); 00267 for (i__ = 1; i__ <= i__2; ++i__) { 00268 a_ref(i__, j) = a_ref(i__, j) * mul; 00269 /* L80: */ 00270 } 00271 /* L90: */ 00272 } 00273 00274 } else if (itype == 4) { 00275 00276 /* Lower half of a symmetric band matrix */ 00277 00278 k3 = *kl + 1; 00279 k4 = *n + 1; 00280 i__1 = *n; 00281 for (j = 1; j <= i__1; ++j) { 00282 /* Computing MIN */ 00283 i__3 = k3, i__4 = k4 - j; 00284 i__2 = minMACRO(i__3,i__4); 00285 for (i__ = 1; i__ <= i__2; ++i__) { 00286 a_ref(i__, j) = a_ref(i__, j) * mul; 00287 /* L100: */ 00288 } 00289 /* L110: */ 00290 } 00291 00292 } else if (itype == 5) { 00293 00294 /* Upper half of a symmetric band matrix */ 00295 00296 k1 = *ku + 2; 00297 k3 = *ku + 1; 00298 i__1 = *n; 00299 for (j = 1; j <= i__1; ++j) { 00300 /* Computing MAX */ 00301 i__2 = k1 - j; 00302 i__3 = k3; 00303 for (i__ = maxMACRO(i__2,1); i__ <= i__3; ++i__) { 00304 a_ref(i__, j) = a_ref(i__, j) * mul; 00305 /* L120: */ 00306 } 00307 /* L130: */ 00308 } 00309 00310 } else if (itype == 6) { 00311 00312 /* Band matrix */ 00313 00314 k1 = *kl + *ku + 2; 00315 k2 = *kl + 1; 00316 k3 = (*kl << 1) + *ku + 1; 00317 k4 = *kl + *ku + 1 + *m; 00318 i__1 = *n; 00319 for (j = 1; j <= i__1; ++j) { 00320 /* Computing MAX */ 00321 i__3 = k1 - j; 00322 /* Computing MIN */ 00323 i__4 = k3, i__5 = k4 - j; 00324 i__2 = minMACRO(i__4,i__5); 00325 for (i__ = maxMACRO(i__3,k2); i__ <= i__2; ++i__) { 00326 a_ref(i__, j) = a_ref(i__, j) * mul; 00327 /* L140: */ 00328 } 00329 /* L150: */ 00330 } 00331 00332 } 00333 00334 if (! done) { 00335 goto L10; 00336 } 00337 00338 return 0; 00339 00340 /* End of DLASCL */ 00341 00342 } /* dlascl_ */ 00343 00344 #undef a_ref 00345 00346 00347 #endif