template_lapack_orgtr.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_ORGTR_HEADER
00036 #define TEMPLATE_LAPACK_ORGTR_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_orgtr(const char *uplo, const integer *n, Treal *a, const integer *
00041         lda, const Treal *tau, Treal *work, const integer *lwork, 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        June 30, 1999   
00047 
00048 
00049     Purpose   
00050     =======   
00051 
00052     DORGTR generates a real orthogonal matrix Q which is defined as the   
00053     product of n-1 elementary reflectors of order N, as returned by   
00054     DSYTRD:   
00055 
00056     if UPLO = 'U', Q = H(n-1) . . . H(2) H(1),   
00057 
00058     if UPLO = 'L', Q = H(1) H(2) . . . H(n-1).   
00059 
00060     Arguments   
00061     =========   
00062 
00063     UPLO    (input) CHARACTER*1   
00064             = 'U': Upper triangle of A contains elementary reflectors   
00065                    from DSYTRD;   
00066             = 'L': Lower triangle of A contains elementary reflectors   
00067                    from DSYTRD.   
00068 
00069     N       (input) INTEGER   
00070             The order of the matrix Q. N >= 0.   
00071 
00072     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)   
00073             On entry, the vectors which define the elementary reflectors,   
00074             as returned by DSYTRD.   
00075             On exit, the N-by-N orthogonal matrix Q.   
00076 
00077     LDA     (input) INTEGER   
00078             The leading dimension of the array A. LDA >= max(1,N).   
00079 
00080     TAU     (input) DOUBLE PRECISION array, dimension (N-1)   
00081             TAU(i) must contain the scalar factor of the elementary   
00082             reflector H(i), as returned by DSYTRD.   
00083 
00084     WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)   
00085             On exit, if INFO = 0, WORK(1) returns the optimal LWORK.   
00086 
00087     LWORK   (input) INTEGER   
00088             The dimension of the array WORK. LWORK >= max(1,N-1).   
00089             For optimum performance LWORK >= (N-1)*NB, where NB is   
00090             the optimal blocksize.   
00091 
00092             If LWORK = -1, then a workspace query is assumed; the routine   
00093             only calculates the optimal size of the WORK array, returns   
00094             this value as the first entry of the WORK array, and no error   
00095             message related to LWORK is issued by XERBLA.   
00096 
00097     INFO    (output) INTEGER   
00098             = 0:  successful exit   
00099             < 0:  if INFO = -i, the i-th argument had an illegal value   
00100 
00101     =====================================================================   
00102 
00103 
00104        Test the input arguments   
00105 
00106        Parameter adjustments */
00107     /* Table of constant values */
00108      integer c__1 = 1;
00109      integer c_n1 = -1;
00110     
00111     /* System generated locals */
00112     integer a_dim1, a_offset, i__1, i__2, i__3;
00113     /* Local variables */
00114      integer i__, j;
00115      integer iinfo;
00116      logical upper;
00117      integer nb;
00118      integer lwkopt;
00119      logical lquery;
00120 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00121 
00122 
00123     a_dim1 = *lda;
00124     a_offset = 1 + a_dim1 * 1;
00125     a -= a_offset;
00126     --tau;
00127     --work;
00128 
00129     /* Initialization added by Elias to get rid of compiler warnings. */
00130     lwkopt = 0;
00131     /* Function Body */
00132     *info = 0;
00133     lquery = *lwork == -1;
00134     upper = template_blas_lsame(uplo, "U");
00135     if (! upper && ! template_blas_lsame(uplo, "L")) {
00136         *info = -1;
00137     } else if (*n < 0) {
00138         *info = -2;
00139     } else if (*lda < maxMACRO(1,*n)) {
00140         *info = -4;
00141     } else /* if(complicated condition) */ {
00142 /* Computing MAX */
00143         i__1 = 1, i__2 = *n - 1;
00144         if (*lwork < maxMACRO(i__1,i__2) && ! lquery) {
00145             *info = -7;
00146         }
00147     }
00148 
00149     if (*info == 0) {
00150         if (upper) {
00151             i__1 = *n - 1;
00152             i__2 = *n - 1;
00153             i__3 = *n - 1;
00154             nb = template_lapack_ilaenv(&c__1, "DORGQL", " ", &i__1, &i__2, &i__3, &c_n1, (
00155                     ftnlen)6, (ftnlen)1);
00156         } else {
00157             i__1 = *n - 1;
00158             i__2 = *n - 1;
00159             i__3 = *n - 1;
00160             nb = template_lapack_ilaenv(&c__1, "DORGQR", " ", &i__1, &i__2, &i__3, &c_n1, (
00161                     ftnlen)6, (ftnlen)1);
00162         }
00163 /* Computing MAX */
00164         i__1 = 1, i__2 = *n - 1;
00165         lwkopt = maxMACRO(i__1,i__2) * nb;
00166         work[1] = (Treal) lwkopt;
00167     }
00168 
00169     if (*info != 0) {
00170         i__1 = -(*info);
00171         template_blas_erbla("ORGTR ", &i__1);
00172         return 0;
00173     } else if (lquery) {
00174         return 0;
00175     }
00176 
00177 /*     Quick return if possible */
00178 
00179     if (*n == 0) {
00180         work[1] = 1.;
00181         return 0;
00182     }
00183 
00184     if (upper) {
00185 
00186 /*        Q was determined by a call to DSYTRD with UPLO = 'U'   
00187 
00188           Shift the vectors which define the elementary reflectors one   
00189           column to the left, and set the last row and column of Q to   
00190           those of the unit matrix */
00191 
00192         i__1 = *n - 1;
00193         for (j = 1; j <= i__1; ++j) {
00194             i__2 = j - 1;
00195             for (i__ = 1; i__ <= i__2; ++i__) {
00196                 a_ref(i__, j) = a_ref(i__, j + 1);
00197 /* L10: */
00198             }
00199             a_ref(*n, j) = 0.;
00200 /* L20: */
00201         }
00202         i__1 = *n - 1;
00203         for (i__ = 1; i__ <= i__1; ++i__) {
00204             a_ref(i__, *n) = 0.;
00205 /* L30: */
00206         }
00207         a_ref(*n, *n) = 1.;
00208 
00209 /*        Generate Q(1:n-1,1:n-1) */
00210 
00211         i__1 = *n - 1;
00212         i__2 = *n - 1;
00213         i__3 = *n - 1;
00214         template_lapack_orgql(&i__1, &i__2, &i__3, &a[a_offset], lda, &tau[1], &work[1], 
00215                 lwork, &iinfo);
00216 
00217     } else {
00218 
00219 /*        Q was determined by a call to DSYTRD with UPLO = 'L'.   
00220 
00221           Shift the vectors which define the elementary reflectors one   
00222           column to the right, and set the first row and column of Q to   
00223           those of the unit matrix */
00224 
00225         for (j = *n; j >= 2; --j) {
00226             a_ref(1, j) = 0.;
00227             i__1 = *n;
00228             for (i__ = j + 1; i__ <= i__1; ++i__) {
00229                 a_ref(i__, j) = a_ref(i__, j - 1);
00230 /* L40: */
00231             }
00232 /* L50: */
00233         }
00234         a_ref(1, 1) = 1.;
00235         i__1 = *n;
00236         for (i__ = 2; i__ <= i__1; ++i__) {
00237             a_ref(i__, 1) = 0.;
00238 /* L60: */
00239         }
00240         if (*n > 1) {
00241 
00242 /*           Generate Q(2:n,2:n) */
00243 
00244             i__1 = *n - 1;
00245             i__2 = *n - 1;
00246             i__3 = *n - 1;
00247             template_lapack_orgqr(
00248                                   &i__1, 
00249                                   &i__2, 
00250                                   &i__3, 
00251                                   &a_ref(2, 2), 
00252                                   lda, 
00253                                   &tau[1], 
00254                                   &work[1],
00255                                   lwork, 
00256                                   &iinfo
00257                                   );
00258         }
00259     }
00260     work[1] = (Treal) lwkopt;
00261     return 0;
00262 
00263 /*     End of DORGTR */
00264 
00265 } /* dorgtr_ */
00266 
00267 #undef a_ref
00268 
00269 
00270 #endif

Generated on Mon Sep 17 14:32:56 2012 for ergo by  doxygen 1.4.7