template_lapack_org2l.h

Go to the documentation of this file.
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

Generated on Wed Nov 21 09:32:01 2012 for ergo by  doxygen 1.4.7