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_ORG2L_HEADER 00036 #define TEMPLATE_LAPACK_ORG2L_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_org2l(const integer *m, const integer *n, const integer *k, Treal * 00041 a, const integer *lda, const Treal *tau, Treal *work, integer *info) 00042 { 00043 /* -- LAPACK routine (version 3.0) -- 00044 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00045 Courant Institute, Argonne National Lab, and Rice University 00046 February 29, 1992 00047 00048 00049 Purpose 00050 ======= 00051 00052 DORG2L generates an m by n real matrix Q with orthonormal columns, 00053 which is defined as the last n columns of a product of k elementary 00054 reflectors of order m 00055 00056 Q = H(k) . . . H(2) H(1) 00057 00058 as returned by DGEQLF. 00059 00060 Arguments 00061 ========= 00062 00063 M (input) INTEGER 00064 The number of rows of the matrix Q. M >= 0. 00065 00066 N (input) INTEGER 00067 The number of columns of the matrix Q. M >= N >= 0. 00068 00069 K (input) INTEGER 00070 The number of elementary reflectors whose product defines the 00071 matrix Q. N >= K >= 0. 00072 00073 A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00074 On entry, the (n-k+i)-th column must contain the vector which 00075 defines the elementary reflector H(i), for i = 1,2,...,k, as 00076 returned by DGEQLF in the last k columns of its array 00077 argument A. 00078 On exit, the m by n matrix Q. 00079 00080 LDA (input) INTEGER 00081 The first dimension of the array A. LDA >= max(1,M). 00082 00083 TAU (input) DOUBLE PRECISION array, dimension (K) 00084 TAU(i) must contain the scalar factor of the elementary 00085 reflector H(i), as returned by DGEQLF. 00086 00087 WORK (workspace) DOUBLE PRECISION array, dimension (N) 00088 00089 INFO (output) INTEGER 00090 = 0: successful exit 00091 < 0: if INFO = -i, the i-th argument has an illegal value 00092 00093 ===================================================================== 00094 00095 00096 Test the input arguments 00097 00098 Parameter adjustments */ 00099 /* Table of constant values */ 00100 integer c__1 = 1; 00101 00102 /* System generated locals */ 00103 integer a_dim1, a_offset, i__1, i__2, i__3; 00104 Treal d__1; 00105 /* Local variables */ 00106 integer i__, j, l; 00107 integer ii; 00108 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00109 00110 00111 a_dim1 = *lda; 00112 a_offset = 1 + a_dim1 * 1; 00113 a -= a_offset; 00114 --tau; 00115 --work; 00116 00117 /* Function Body */ 00118 *info = 0; 00119 if (*m < 0) { 00120 *info = -1; 00121 } else if (*n < 0 || *n > *m) { 00122 *info = -2; 00123 } else if (*k < 0 || *k > *n) { 00124 *info = -3; 00125 } else if (*lda < maxMACRO(1,*m)) { 00126 *info = -5; 00127 } 00128 if (*info != 0) { 00129 i__1 = -(*info); 00130 template_blas_erbla("ORG2L ", &i__1); 00131 return 0; 00132 } 00133 00134 /* Quick return if possible */ 00135 00136 if (*n <= 0) { 00137 return 0; 00138 } 00139 00140 /* Initialise columns 1:n-k to columns of the unit matrix */ 00141 00142 i__1 = *n - *k; 00143 for (j = 1; j <= i__1; ++j) { 00144 i__2 = *m; 00145 for (l = 1; l <= i__2; ++l) { 00146 a_ref(l, j) = 0.; 00147 /* L10: */ 00148 } 00149 a_ref(*m - *n + j, j) = 1.; 00150 /* L20: */ 00151 } 00152 00153 i__1 = *k; 00154 for (i__ = 1; i__ <= i__1; ++i__) { 00155 ii = *n - *k + i__; 00156 00157 /* Apply H(i) to A(1:m-k+i,1:n-k+i) from the left */ 00158 00159 a_ref(*m - *n + ii, ii) = 1.; 00160 i__2 = *m - *n + ii; 00161 i__3 = ii - 1; 00162 template_lapack_larf("Left", &i__2, &i__3, &a_ref(1, ii), &c__1, &tau[i__], &a[ 00163 a_offset], lda, &work[1]); 00164 i__2 = *m - *n + ii - 1; 00165 d__1 = -tau[i__]; 00166 template_blas_scal(&i__2, &d__1, &a_ref(1, ii), &c__1); 00167 a_ref(*m - *n + ii, ii) = 1. - tau[i__]; 00168 00169 /* Set A(m-k+i+1:m,n-k+i) to zero */ 00170 00171 i__2 = *m; 00172 for (l = *m - *n + ii + 1; l <= i__2; ++l) { 00173 a_ref(l, ii) = 0.; 00174 /* L30: */ 00175 } 00176 /* L40: */ 00177 } 00178 return 0; 00179 00180 /* End of DORG2L */ 00181 00182 } /* dorg2l_ */ 00183 00184 #undef a_ref 00185 00186 00187 #endif