template_lapack_sytd2.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_SYTD2_HEADER
00036 #define TEMPLATE_LAPACK_SYTD2_HEADER
00037 
00038 #include "template_lapack_common.h"
00039 
00040 template<class Treal>
00041 int template_lapack_sytd2(const char *uplo, const integer *n, Treal *a, const integer *
00042         lda, Treal *d__, Treal *e, Treal *tau, 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        October 31, 1992   
00048 
00049 
00050     Purpose   
00051     =======   
00052 
00053     DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal   
00054     form T by an orthogonal similarity transformation: Q' * A * Q = T.   
00055 
00056     Arguments   
00057     =========   
00058 
00059     UPLO    (input) CHARACTER*1   
00060             Specifies whether the upper or lower triangular part of the   
00061             symmetric matrix A is stored:   
00062             = 'U':  Upper triangular   
00063             = 'L':  Lower triangular   
00064 
00065     N       (input) INTEGER   
00066             The order of the matrix A.  N >= 0.   
00067 
00068     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)   
00069             On entry, the symmetric matrix A.  If UPLO = 'U', the leading   
00070             n-by-n upper triangular part of A contains the upper   
00071             triangular part of the matrix A, and the strictly lower   
00072             triangular part of A is not referenced.  If UPLO = 'L', the   
00073             leading n-by-n lower triangular part of A contains the lower   
00074             triangular part of the matrix A, and the strictly upper   
00075             triangular part of A is not referenced.   
00076             On exit, if UPLO = 'U', the diagonal and first superdiagonal   
00077             of A are overwritten by the corresponding elements of the   
00078             tridiagonal matrix T, and the elements above the first   
00079             superdiagonal, with the array TAU, represent the orthogonal   
00080             matrix Q as a product of elementary reflectors; if UPLO   
00081             = 'L', the diagonal and first subdiagonal of A are over-   
00082             written by the corresponding elements of the tridiagonal   
00083             matrix T, and the elements below the first subdiagonal, with   
00084             the array TAU, represent the orthogonal matrix Q as a product   
00085             of elementary reflectors. See Further Details.   
00086 
00087     LDA     (input) INTEGER   
00088             The leading dimension of the array A.  LDA >= max(1,N).   
00089 
00090     D       (output) DOUBLE PRECISION array, dimension (N)   
00091             The diagonal elements of the tridiagonal matrix T:   
00092             D(i) = A(i,i).   
00093 
00094     E       (output) DOUBLE PRECISION array, dimension (N-1)   
00095             The off-diagonal elements of the tridiagonal matrix T:   
00096             E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.   
00097 
00098     TAU     (output) DOUBLE PRECISION array, dimension (N-1)   
00099             The scalar factors of the elementary reflectors (see Further   
00100             Details).   
00101 
00102     INFO    (output) INTEGER   
00103             = 0:  successful exit   
00104             < 0:  if INFO = -i, the i-th argument had an illegal value.   
00105 
00106     Further Details   
00107     ===============   
00108 
00109     If UPLO = 'U', the matrix Q is represented as a product of elementary   
00110     reflectors   
00111 
00112        Q = H(n-1) . . . H(2) H(1).   
00113 
00114     Each H(i) has the form   
00115 
00116        H(i) = I - tau * v * v'   
00117 
00118     where tau is a real scalar, and v is a real vector with   
00119     v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in   
00120     A(1:i-1,i+1), and tau in TAU(i).   
00121 
00122     If UPLO = 'L', the matrix Q is represented as a product of elementary   
00123     reflectors   
00124 
00125        Q = H(1) H(2) . . . H(n-1).   
00126 
00127     Each H(i) has the form   
00128 
00129        H(i) = I - tau * v * v'   
00130 
00131     where tau is a real scalar, and v is a real vector with   
00132     v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),   
00133     and tau in TAU(i).   
00134 
00135     The contents of A on exit are illustrated by the following examples   
00136     with n = 5:   
00137 
00138     if UPLO = 'U':                       if UPLO = 'L':   
00139 
00140       (  d   e   v2  v3  v4 )              (  d                  )   
00141       (      d   e   v3  v4 )              (  e   d              )   
00142       (          d   e   v4 )              (  v1  e   d          )   
00143       (              d   e  )              (  v1  v2  e   d      )   
00144       (                  d  )              (  v1  v2  v3  e   d  )   
00145 
00146     where d and e denote diagonal and off-diagonal elements of T, and vi   
00147     denotes an element of the vector defining H(i).   
00148 
00149     =====================================================================   
00150 
00151 
00152        Test the input parameters   
00153 
00154        Parameter adjustments */
00155     /* Table of constant values */
00156      integer c__1 = 1;
00157      Treal c_b8 = 0.;
00158      Treal c_b14 = -1.;
00159     
00160     /* System generated locals */
00161     integer a_dim1, a_offset, i__1, i__2, i__3;
00162     /* Local variables */
00163      Treal taui;
00164      integer i__;
00165      Treal alpha;
00166      logical upper;
00167 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00168 
00169 
00170     a_dim1 = *lda;
00171     a_offset = 1 + a_dim1 * 1;
00172     a -= a_offset;
00173     --d__;
00174     --e;
00175     --tau;
00176 
00177     /* Function Body */
00178     *info = 0;
00179     upper = template_blas_lsame(uplo, "U");
00180     if (! upper && ! template_blas_lsame(uplo, "L")) {
00181         *info = -1;
00182     } else if (*n < 0) {
00183         *info = -2;
00184     } else if (*lda < maxMACRO(1,*n)) {
00185         *info = -4;
00186     }
00187     if (*info != 0) {
00188         i__1 = -(*info);
00189         template_blas_erbla("SYTD2 ", &i__1);
00190         return 0;
00191     }
00192 
00193 /*     Quick return if possible */
00194 
00195     if (*n <= 0) {
00196         return 0;
00197     }
00198 
00199     if (upper) {
00200 
00201 /*        Reduce the upper triangle of A */
00202 
00203         for (i__ = *n - 1; i__ >= 1; --i__) {
00204 
00205 /*           Generate elementary reflector H(i) = I - tau * v * v'   
00206              to annihilate A(1:i-1,i+1) */
00207 
00208             template_lapack_larfg(&i__, &a_ref(i__, i__ + 1), &a_ref(1, i__ + 1), &c__1, &
00209                     taui);
00210             e[i__] = a_ref(i__, i__ + 1);
00211 
00212             if (taui != 0.) {
00213 
00214 /*              Apply H(i) from both sides to A(1:i,1:i) */
00215 
00216                 a_ref(i__, i__ + 1) = 1.;
00217 
00218 /*              Compute  x := tau * A * v  storing x in TAU(1:i) */
00219 
00220                 template_blas_symv(uplo, &i__, &taui, &a[a_offset], lda, &a_ref(1, i__ + 
00221                         1), &c__1, &c_b8, &tau[1], &c__1);
00222 
00223 /*              Compute  w := x - 1/2 * tau * (x'*v) * v */
00224 
00225                 alpha = taui * -.5 * template_blas_dot(&i__, &tau[1], &c__1, &a_ref(1, 
00226                         i__ + 1), &c__1);
00227                 template_blas_axpy(&i__, &alpha, &a_ref(1, i__ + 1), &c__1, &tau[1], &
00228                         c__1);
00229 
00230 /*              Apply the transformation as a rank-2 update:   
00231                    A := A - v * w' - w * v' */
00232 
00233                 template_blas_syr2(uplo, &i__, &c_b14, &a_ref(1, i__ + 1), &c__1, &tau[1],
00234                          &c__1, &a[a_offset], lda);
00235 
00236                 a_ref(i__, i__ + 1) = e[i__];
00237             }
00238             d__[i__ + 1] = a_ref(i__ + 1, i__ + 1);
00239             tau[i__] = taui;
00240 /* L10: */
00241         }
00242         d__[1] = a_ref(1, 1);
00243     } else {
00244 
00245 /*        Reduce the lower triangle of A */
00246 
00247         i__1 = *n - 1;
00248         for (i__ = 1; i__ <= i__1; ++i__) {
00249 
00250 /*           Generate elementary reflector H(i) = I - tau * v * v'   
00251              to annihilate A(i+2:n,i)   
00252 
00253    Computing MIN */
00254             i__2 = i__ + 2;
00255             i__3 = *n - i__;
00256             template_lapack_larfg(&i__3, &a_ref(i__ + 1, i__), &a_ref(minMACRO(i__2,*n), i__), &
00257                     c__1, &taui);
00258             e[i__] = a_ref(i__ + 1, i__);
00259 
00260             if (taui != 0.) {
00261 
00262 /*              Apply H(i) from both sides to A(i+1:n,i+1:n) */
00263 
00264                 a_ref(i__ + 1, i__) = 1.;
00265 
00266 /*              Compute  x := tau * A * v  storing y in TAU(i:n-1) */
00267 
00268                 i__2 = *n - i__;
00269                 template_blas_symv(uplo, &i__2, &taui, &a_ref(i__ + 1, i__ + 1), lda, &
00270                         a_ref(i__ + 1, i__), &c__1, &c_b8, &tau[i__], &c__1);
00271 
00272 /*              Compute  w := x - 1/2 * tau * (x'*v) * v */
00273 
00274                 i__2 = *n - i__;
00275                 alpha = taui * -.5 * template_blas_dot(&i__2, &tau[i__], &c__1, &a_ref(
00276                         i__ + 1, i__), &c__1);
00277                 i__2 = *n - i__;
00278                 template_blas_axpy(&i__2, &alpha, &a_ref(i__ + 1, i__), &c__1, &tau[i__], 
00279                         &c__1);
00280 
00281 /*              Apply the transformation as a rank-2 update:   
00282                    A := A - v * w' - w * v' */
00283 
00284                 i__2 = *n - i__;
00285                 template_blas_syr2(uplo, &i__2, &c_b14, &a_ref(i__ + 1, i__), &c__1, &tau[
00286                         i__], &c__1, &a_ref(i__ + 1, i__ + 1), lda)
00287                         ;
00288 
00289                 a_ref(i__ + 1, i__) = e[i__];
00290             }
00291             d__[i__] = a_ref(i__, i__);
00292             tau[i__] = taui;
00293 /* L20: */
00294         }
00295         d__[*n] = a_ref(*n, *n);
00296     }
00297 
00298     return 0;
00299 
00300 /*     End of DSYTD2 */
00301 
00302 } /* dsytd2_ */
00303 
00304 #undef a_ref
00305 
00306 
00307 #endif

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1