template_lapack_orm2r.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_ORM2R_HEADER
00036 #define TEMPLATE_LAPACK_ORM2R_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_orm2r(const char *side, const char *trans, const integer *m, const integer *n, 
00041         const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal *
00042         c__, const integer *ldc, Treal *work, integer *info)
00043 {
00044 /*  -- LAPACK 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     DORM2R overwrites the general real m by n matrix C with   
00054 
00055           Q * C  if SIDE = 'L' and TRANS = 'N', or   
00056 
00057           Q'* C  if SIDE = 'L' and TRANS = 'T', or   
00058 
00059           C * Q  if SIDE = 'R' and TRANS = 'N', or   
00060 
00061           C * Q' if SIDE = 'R' and TRANS = 'T',   
00062 
00063     where Q is a real orthogonal matrix defined as the product of k   
00064     elementary reflectors   
00065 
00066           Q = H(1) H(2) . . . H(k)   
00067 
00068     as returned by DGEQRF. Q is of order m if SIDE = 'L' and of order n   
00069     if SIDE = 'R'.   
00070 
00071     Arguments   
00072     =========   
00073 
00074     SIDE    (input) CHARACTER*1   
00075             = 'L': apply Q or Q' from the Left   
00076             = 'R': apply Q or Q' from the Right   
00077 
00078     TRANS   (input) CHARACTER*1   
00079             = 'N': apply Q  (No transpose)   
00080             = 'T': apply Q' (Transpose)   
00081 
00082     M       (input) INTEGER   
00083             The number of rows of the matrix C. M >= 0.   
00084 
00085     N       (input) INTEGER   
00086             The number of columns of the matrix C. N >= 0.   
00087 
00088     K       (input) INTEGER   
00089             The number of elementary reflectors whose product defines   
00090             the matrix Q.   
00091             If SIDE = 'L', M >= K >= 0;   
00092             if SIDE = 'R', N >= K >= 0.   
00093 
00094     A       (input) DOUBLE PRECISION array, dimension (LDA,K)   
00095             The i-th column must contain the vector which defines the   
00096             elementary reflector H(i), for i = 1,2,...,k, as returned by   
00097             DGEQRF in the first k columns of its array argument A.   
00098             A is modified by the routine but restored on exit.   
00099 
00100     LDA     (input) INTEGER   
00101             The leading dimension of the array A.   
00102             If SIDE = 'L', LDA >= max(1,M);   
00103             if SIDE = 'R', LDA >= max(1,N).   
00104 
00105     TAU     (input) DOUBLE PRECISION array, dimension (K)   
00106             TAU(i) must contain the scalar factor of the elementary   
00107             reflector H(i), as returned by DGEQRF.   
00108 
00109     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)   
00110             On entry, the m by n matrix C.   
00111             On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q.   
00112 
00113     LDC     (input) INTEGER   
00114             The leading dimension of the array C. LDC >= max(1,M).   
00115 
00116     WORK    (workspace) DOUBLE PRECISION array, dimension   
00117                                      (N) if SIDE = 'L',   
00118                                      (M) if SIDE = 'R'   
00119 
00120     INFO    (output) INTEGER   
00121             = 0: successful exit   
00122             < 0: if INFO = -i, the i-th argument had an illegal value   
00123 
00124     =====================================================================   
00125 
00126 
00127        Test the input arguments   
00128 
00129        Parameter adjustments */
00130     /* Table of constant values */
00131      integer c__1 = 1;
00132     
00133     /* System generated locals */
00134     integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2;
00135     /* Local variables */
00136      logical left;
00137      integer i__;
00138      integer i1, i2, i3, ic, jc, mi, ni, nq;
00139      logical notran;
00140      Treal aii;
00141 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00142 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
00143 
00144 
00145     a_dim1 = *lda;
00146     a_offset = 1 + a_dim1 * 1;
00147     a -= a_offset;
00148     --tau;
00149     c_dim1 = *ldc;
00150     c_offset = 1 + c_dim1 * 1;
00151     c__ -= c_offset;
00152     --work;
00153 
00154     /* Function Body */
00155     *info = 0;
00156     left = template_blas_lsame(side, "L");
00157     notran = template_blas_lsame(trans, "N");
00158 
00159 /*     NQ is the order of Q */
00160 
00161     if (left) {
00162         nq = *m;
00163     } else {
00164         nq = *n;
00165     }
00166     if (! left && ! template_blas_lsame(side, "R")) {
00167         *info = -1;
00168     } else if (! notran && ! template_blas_lsame(trans, "T")) {
00169         *info = -2;
00170     } else if (*m < 0) {
00171         *info = -3;
00172     } else if (*n < 0) {
00173         *info = -4;
00174     } else if (*k < 0 || *k > nq) {
00175         *info = -5;
00176     } else if (*lda < maxMACRO(1,nq)) {
00177         *info = -7;
00178     } else if (*ldc < maxMACRO(1,*m)) {
00179         *info = -10;
00180     }
00181     if (*info != 0) {
00182         i__1 = -(*info);
00183         template_blas_erbla("ORM2R ", &i__1);
00184         return 0;
00185     }
00186 
00187 /*     Quick return if possible */
00188 
00189     if (*m == 0 || *n == 0 || *k == 0) {
00190         return 0;
00191     }
00192 
00193     if ( ( left && ! notran ) || ( ! left && notran ) ) {
00194         i1 = 1;
00195         i2 = *k;
00196         i3 = 1;
00197     } else {
00198         i1 = *k;
00199         i2 = 1;
00200         i3 = -1;
00201     }
00202 
00203     if (left) {
00204         ni = *n;
00205         jc = 1;
00206     } else {
00207         mi = *m;
00208         ic = 1;
00209     }
00210 
00211     i__1 = i2;
00212     i__2 = i3;
00213     for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
00214         if (left) {
00215 
00216 /*           H(i) is applied to C(i:m,1:n) */
00217 
00218             mi = *m - i__ + 1;
00219             ic = i__;
00220         } else {
00221 
00222 /*           H(i) is applied to C(1:m,i:n) */
00223 
00224             ni = *n - i__ + 1;
00225             jc = i__;
00226         }
00227 
00228 /*        Apply H(i) */
00229 
00230         aii = a_ref(i__, i__);
00231         a_ref(i__, i__) = 1.;
00232         template_lapack_larf(side, &mi, &ni, &a_ref(i__, i__), &c__1, &tau[i__], &c___ref(
00233                 ic, jc), ldc, &work[1]);
00234         a_ref(i__, i__) = aii;
00235 /* L10: */
00236     }
00237     return 0;
00238 
00239 /*     End of DORM2R */
00240 
00241 } /* dorm2r_ */
00242 
00243 #undef c___ref
00244 #undef a_ref
00245 
00246 
00247 #endif

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