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_LASR_HEADER 00036 #define TEMPLATE_LAPACK_LASR_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_lasr(const char *side, const char *pivot, const char *direct, const integer *m, 00041 const integer *n, const Treal *c__, const Treal *s, Treal *a, const integer * 00042 lda) 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 DLASR performs the transformation 00054 00055 A := P*A, when SIDE = 'L' or 'l' ( Left-hand side ) 00056 00057 A := A*P', when SIDE = 'R' or 'r' ( Right-hand side ) 00058 00059 where A is an m by n real matrix and P is an orthogonal matrix, 00060 consisting of a sequence of plane rotations determined by the 00061 parameters PIVOT and DIRECT as follows ( z = m when SIDE = 'L' or 'l' 00062 and z = n when SIDE = 'R' or 'r' ): 00063 00064 When DIRECT = 'F' or 'f' ( Forward sequence ) then 00065 00066 P = P( z - 1 )*...*P( 2 )*P( 1 ), 00067 00068 and when DIRECT = 'B' or 'b' ( Backward sequence ) then 00069 00070 P = P( 1 )*P( 2 )*...*P( z - 1 ), 00071 00072 where P( k ) is a plane rotation matrix for the following planes: 00073 00074 when PIVOT = 'V' or 'v' ( Variable pivot ), 00075 the plane ( k, k + 1 ) 00076 00077 when PIVOT = 'T' or 't' ( Top pivot ), 00078 the plane ( 1, k + 1 ) 00079 00080 when PIVOT = 'B' or 'b' ( Bottom pivot ), 00081 the plane ( k, z ) 00082 00083 c( k ) and s( k ) must contain the cosine and sine that define the 00084 matrix P( k ). The two by two plane rotation part of the matrix 00085 P( k ), R( k ), is assumed to be of the form 00086 00087 R( k ) = ( c( k ) s( k ) ). 00088 ( -s( k ) c( k ) ) 00089 00090 This version vectorises across rows of the array A when SIDE = 'L'. 00091 00092 Arguments 00093 ========= 00094 00095 SIDE (input) CHARACTER*1 00096 Specifies whether the plane rotation matrix P is applied to 00097 A on the left or the right. 00098 = 'L': Left, compute A := P*A 00099 = 'R': Right, compute A:= A*P' 00100 00101 DIRECT (input) CHARACTER*1 00102 Specifies whether P is a forward or backward sequence of 00103 plane rotations. 00104 = 'F': Forward, P = P( z - 1 )*...*P( 2 )*P( 1 ) 00105 = 'B': Backward, P = P( 1 )*P( 2 )*...*P( z - 1 ) 00106 00107 PIVOT (input) CHARACTER*1 00108 Specifies the plane for which P(k) is a plane rotation 00109 matrix. 00110 = 'V': Variable pivot, the plane (k,k+1) 00111 = 'T': Top pivot, the plane (1,k+1) 00112 = 'B': Bottom pivot, the plane (k,z) 00113 00114 M (input) INTEGER 00115 The number of rows of the matrix A. If m <= 1, an immediate 00116 return is effected. 00117 00118 N (input) INTEGER 00119 The number of columns of the matrix A. If n <= 1, an 00120 immediate return is effected. 00121 00122 C, S (input) DOUBLE PRECISION arrays, dimension 00123 (M-1) if SIDE = 'L' 00124 (N-1) if SIDE = 'R' 00125 c(k) and s(k) contain the cosine and sine that define the 00126 matrix P(k). The two by two plane rotation part of the 00127 matrix P(k), R(k), is assumed to be of the form 00128 R( k ) = ( c( k ) s( k ) ). 00129 ( -s( k ) c( k ) ) 00130 00131 A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00132 The m by n matrix A. On exit, A is overwritten by P*A if 00133 SIDE = 'R' or by A*P' if SIDE = 'L'. 00134 00135 LDA (input) INTEGER 00136 The leading dimension of the array A. LDA >= max(1,M). 00137 00138 ===================================================================== 00139 00140 00141 Test the input parameters 00142 00143 Parameter adjustments */ 00144 /* System generated locals */ 00145 integer a_dim1, a_offset, i__1, i__2; 00146 /* Local variables */ 00147 integer info; 00148 Treal temp; 00149 integer i__, j; 00150 Treal ctemp, stemp; 00151 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00152 00153 --c__; 00154 --s; 00155 a_dim1 = *lda; 00156 a_offset = 1 + a_dim1 * 1; 00157 a -= a_offset; 00158 00159 /* Function Body */ 00160 info = 0; 00161 if (! (template_blas_lsame(side, "L") || template_blas_lsame(side, "R"))) { 00162 info = 1; 00163 } else if (! (template_blas_lsame(pivot, "V") || template_blas_lsame(pivot, 00164 "T") || template_blas_lsame(pivot, "B"))) { 00165 info = 2; 00166 } else if (! (template_blas_lsame(direct, "F") || template_blas_lsame(direct, 00167 "B"))) { 00168 info = 3; 00169 } else if (*m < 0) { 00170 info = 4; 00171 } else if (*n < 0) { 00172 info = 5; 00173 } else if (*lda < maxMACRO(1,*m)) { 00174 info = 9; 00175 } 00176 if (info != 0) { 00177 template_blas_erbla("LASR ", &info); 00178 return 0; 00179 } 00180 00181 /* Quick return if possible */ 00182 00183 if (*m == 0 || *n == 0) { 00184 return 0; 00185 } 00186 if (template_blas_lsame(side, "L")) { 00187 00188 /* Form P * A */ 00189 00190 if (template_blas_lsame(pivot, "V")) { 00191 if (template_blas_lsame(direct, "F")) { 00192 i__1 = *m - 1; 00193 for (j = 1; j <= i__1; ++j) { 00194 ctemp = c__[j]; 00195 stemp = s[j]; 00196 if (ctemp != 1. || stemp != 0.) { 00197 i__2 = *n; 00198 for (i__ = 1; i__ <= i__2; ++i__) { 00199 temp = a_ref(j + 1, i__); 00200 a_ref(j + 1, i__) = ctemp * temp - stemp * a_ref( 00201 j, i__); 00202 a_ref(j, i__) = stemp * temp + ctemp * a_ref(j, 00203 i__); 00204 /* L10: */ 00205 } 00206 } 00207 /* L20: */ 00208 } 00209 } else if (template_blas_lsame(direct, "B")) { 00210 for (j = *m - 1; j >= 1; --j) { 00211 ctemp = c__[j]; 00212 stemp = s[j]; 00213 if (ctemp != 1. || stemp != 0.) { 00214 i__1 = *n; 00215 for (i__ = 1; i__ <= i__1; ++i__) { 00216 temp = a_ref(j + 1, i__); 00217 a_ref(j + 1, i__) = ctemp * temp - stemp * a_ref( 00218 j, i__); 00219 a_ref(j, i__) = stemp * temp + ctemp * a_ref(j, 00220 i__); 00221 /* L30: */ 00222 } 00223 } 00224 /* L40: */ 00225 } 00226 } 00227 } else if (template_blas_lsame(pivot, "T")) { 00228 if (template_blas_lsame(direct, "F")) { 00229 i__1 = *m; 00230 for (j = 2; j <= i__1; ++j) { 00231 ctemp = c__[j - 1]; 00232 stemp = s[j - 1]; 00233 if (ctemp != 1. || stemp != 0.) { 00234 i__2 = *n; 00235 for (i__ = 1; i__ <= i__2; ++i__) { 00236 temp = a_ref(j, i__); 00237 a_ref(j, i__) = ctemp * temp - stemp * a_ref(1, 00238 i__); 00239 a_ref(1, i__) = stemp * temp + ctemp * a_ref(1, 00240 i__); 00241 /* L50: */ 00242 } 00243 } 00244 /* L60: */ 00245 } 00246 } else if (template_blas_lsame(direct, "B")) { 00247 for (j = *m; j >= 2; --j) { 00248 ctemp = c__[j - 1]; 00249 stemp = s[j - 1]; 00250 if (ctemp != 1. || stemp != 0.) { 00251 i__1 = *n; 00252 for (i__ = 1; i__ <= i__1; ++i__) { 00253 temp = a_ref(j, i__); 00254 a_ref(j, i__) = ctemp * temp - stemp * a_ref(1, 00255 i__); 00256 a_ref(1, i__) = stemp * temp + ctemp * a_ref(1, 00257 i__); 00258 /* L70: */ 00259 } 00260 } 00261 /* L80: */ 00262 } 00263 } 00264 } else if (template_blas_lsame(pivot, "B")) { 00265 if (template_blas_lsame(direct, "F")) { 00266 i__1 = *m - 1; 00267 for (j = 1; j <= i__1; ++j) { 00268 ctemp = c__[j]; 00269 stemp = s[j]; 00270 if (ctemp != 1. || stemp != 0.) { 00271 i__2 = *n; 00272 for (i__ = 1; i__ <= i__2; ++i__) { 00273 temp = a_ref(j, i__); 00274 a_ref(j, i__) = stemp * a_ref(*m, i__) + ctemp * 00275 temp; 00276 a_ref(*m, i__) = ctemp * a_ref(*m, i__) - stemp * 00277 temp; 00278 /* L90: */ 00279 } 00280 } 00281 /* L100: */ 00282 } 00283 } else if (template_blas_lsame(direct, "B")) { 00284 for (j = *m - 1; j >= 1; --j) { 00285 ctemp = c__[j]; 00286 stemp = s[j]; 00287 if (ctemp != 1. || stemp != 0.) { 00288 i__1 = *n; 00289 for (i__ = 1; i__ <= i__1; ++i__) { 00290 temp = a_ref(j, i__); 00291 a_ref(j, i__) = stemp * a_ref(*m, i__) + ctemp * 00292 temp; 00293 a_ref(*m, i__) = ctemp * a_ref(*m, i__) - stemp * 00294 temp; 00295 /* L110: */ 00296 } 00297 } 00298 /* L120: */ 00299 } 00300 } 00301 } 00302 } else if (template_blas_lsame(side, "R")) { 00303 00304 /* Form A * P' */ 00305 00306 if (template_blas_lsame(pivot, "V")) { 00307 if (template_blas_lsame(direct, "F")) { 00308 i__1 = *n - 1; 00309 for (j = 1; j <= i__1; ++j) { 00310 ctemp = c__[j]; 00311 stemp = s[j]; 00312 if (ctemp != 1. || stemp != 0.) { 00313 i__2 = *m; 00314 for (i__ = 1; i__ <= i__2; ++i__) { 00315 temp = a_ref(i__, j + 1); 00316 a_ref(i__, j + 1) = ctemp * temp - stemp * a_ref( 00317 i__, j); 00318 a_ref(i__, j) = stemp * temp + ctemp * a_ref(i__, 00319 j); 00320 /* L130: */ 00321 } 00322 } 00323 /* L140: */ 00324 } 00325 } else if (template_blas_lsame(direct, "B")) { 00326 for (j = *n - 1; j >= 1; --j) { 00327 ctemp = c__[j]; 00328 stemp = s[j]; 00329 if (ctemp != 1. || stemp != 0.) { 00330 i__1 = *m; 00331 for (i__ = 1; i__ <= i__1; ++i__) { 00332 temp = a_ref(i__, j + 1); 00333 a_ref(i__, j + 1) = ctemp * temp - stemp * a_ref( 00334 i__, j); 00335 a_ref(i__, j) = stemp * temp + ctemp * a_ref(i__, 00336 j); 00337 /* L150: */ 00338 } 00339 } 00340 /* L160: */ 00341 } 00342 } 00343 } else if (template_blas_lsame(pivot, "T")) { 00344 if (template_blas_lsame(direct, "F")) { 00345 i__1 = *n; 00346 for (j = 2; j <= i__1; ++j) { 00347 ctemp = c__[j - 1]; 00348 stemp = s[j - 1]; 00349 if (ctemp != 1. || stemp != 0.) { 00350 i__2 = *m; 00351 for (i__ = 1; i__ <= i__2; ++i__) { 00352 temp = a_ref(i__, j); 00353 a_ref(i__, j) = ctemp * temp - stemp * a_ref(i__, 00354 1); 00355 a_ref(i__, 1) = stemp * temp + ctemp * a_ref(i__, 00356 1); 00357 /* L170: */ 00358 } 00359 } 00360 /* L180: */ 00361 } 00362 } else if (template_blas_lsame(direct, "B")) { 00363 for (j = *n; j >= 2; --j) { 00364 ctemp = c__[j - 1]; 00365 stemp = s[j - 1]; 00366 if (ctemp != 1. || stemp != 0.) { 00367 i__1 = *m; 00368 for (i__ = 1; i__ <= i__1; ++i__) { 00369 temp = a_ref(i__, j); 00370 a_ref(i__, j) = ctemp * temp - stemp * a_ref(i__, 00371 1); 00372 a_ref(i__, 1) = stemp * temp + ctemp * a_ref(i__, 00373 1); 00374 /* L190: */ 00375 } 00376 } 00377 /* L200: */ 00378 } 00379 } 00380 } else if (template_blas_lsame(pivot, "B")) { 00381 if (template_blas_lsame(direct, "F")) { 00382 i__1 = *n - 1; 00383 for (j = 1; j <= i__1; ++j) { 00384 ctemp = c__[j]; 00385 stemp = s[j]; 00386 if (ctemp != 1. || stemp != 0.) { 00387 i__2 = *m; 00388 for (i__ = 1; i__ <= i__2; ++i__) { 00389 temp = a_ref(i__, j); 00390 a_ref(i__, j) = stemp * a_ref(i__, *n) + ctemp * 00391 temp; 00392 a_ref(i__, *n) = ctemp * a_ref(i__, *n) - stemp * 00393 temp; 00394 /* L210: */ 00395 } 00396 } 00397 /* L220: */ 00398 } 00399 } else if (template_blas_lsame(direct, "B")) { 00400 for (j = *n - 1; j >= 1; --j) { 00401 ctemp = c__[j]; 00402 stemp = s[j]; 00403 if (ctemp != 1. || stemp != 0.) { 00404 i__1 = *m; 00405 for (i__ = 1; i__ <= i__1; ++i__) { 00406 temp = a_ref(i__, j); 00407 a_ref(i__, j) = stemp * a_ref(i__, *n) + ctemp * 00408 temp; 00409 a_ref(i__, *n) = ctemp * a_ref(i__, *n) - stemp * 00410 temp; 00411 /* L230: */ 00412 } 00413 } 00414 /* L240: */ 00415 } 00416 } 00417 } 00418 } 00419 00420 return 0; 00421 00422 /* End of DLASR */ 00423 00424 } /* dlasr_ */ 00425 00426 #undef a_ref 00427 00428 00429 #endif