template_lapack_larft.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_LARFT_HEADER
00036 #define TEMPLATE_LAPACK_LARFT_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_larft(const char *direct, const char *storev, const integer *n, const integer *
00041         k, Treal *v, const integer *ldv, const Treal *tau, Treal *t, 
00042         const integer *ldt)
00043 {
00044 /*  -- LAPACK auxiliary 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     DLARFT forms the triangular factor T of a real block reflector H   
00054     of order n, which is defined as a product of k elementary reflectors.   
00055 
00056     If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;   
00057 
00058     If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.   
00059 
00060     If STOREV = 'C', the vector which defines the elementary reflector   
00061     H(i) is stored in the i-th column of the array V, and   
00062 
00063        H  =  I - V * T * V'   
00064 
00065     If STOREV = 'R', the vector which defines the elementary reflector   
00066     H(i) is stored in the i-th row of the array V, and   
00067 
00068        H  =  I - V' * T * V   
00069 
00070     Arguments   
00071     =========   
00072 
00073     DIRECT  (input) CHARACTER*1   
00074             Specifies the order in which the elementary reflectors are   
00075             multiplied to form the block reflector:   
00076             = 'F': H = H(1) H(2) . . . H(k) (Forward)   
00077             = 'B': H = H(k) . . . H(2) H(1) (Backward)   
00078 
00079     STOREV  (input) CHARACTER*1   
00080             Specifies how the vectors which define the elementary   
00081             reflectors are stored (see also Further Details):   
00082             = 'C': columnwise   
00083             = 'R': rowwise   
00084 
00085     N       (input) INTEGER   
00086             The order of the block reflector H. N >= 0.   
00087 
00088     K       (input) INTEGER   
00089             The order of the triangular factor T (= the number of   
00090             elementary reflectors). K >= 1.   
00091 
00092     V       (input/output) DOUBLE PRECISION array, dimension   
00093                                  (LDV,K) if STOREV = 'C'   
00094                                  (LDV,N) if STOREV = 'R'   
00095             The matrix V. See further details.   
00096 
00097     LDV     (input) INTEGER   
00098             The leading dimension of the array V.   
00099             If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.   
00100 
00101     TAU     (input) DOUBLE PRECISION array, dimension (K)   
00102             TAU(i) must contain the scalar factor of the elementary   
00103             reflector H(i).   
00104 
00105     T       (output) DOUBLE PRECISION array, dimension (LDT,K)   
00106             The k by k triangular factor T of the block reflector.   
00107             If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is   
00108             lower triangular. The rest of the array is not used.   
00109 
00110     LDT     (input) INTEGER   
00111             The leading dimension of the array T. LDT >= K.   
00112 
00113     Further Details   
00114     ===============   
00115 
00116     The shape of the matrix V and the storage of the vectors which define   
00117     the H(i) is best illustrated by the following example with n = 5 and   
00118     k = 3. The elements equal to 1 are not stored; the corresponding   
00119     array elements are modified but restored on exit. The rest of the   
00120     array is not used.   
00121 
00122     DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':   
00123 
00124                  V = (  1       )                 V = (  1 v1 v1 v1 v1 )   
00125                      ( v1  1    )                     (     1 v2 v2 v2 )   
00126                      ( v1 v2  1 )                     (        1 v3 v3 )   
00127                      ( v1 v2 v3 )   
00128                      ( v1 v2 v3 )   
00129 
00130     DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':   
00131 
00132                  V = ( v1 v2 v3 )                 V = ( v1 v1  1       )   
00133                      ( v1 v2 v3 )                     ( v2 v2 v2  1    )   
00134                      (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )   
00135                      (     1 v3 )   
00136                      (        1 )   
00137 
00138     =====================================================================   
00139 
00140 
00141        Quick return if possible   
00142 
00143        Parameter adjustments */
00144     /* Table of constant values */
00145      integer c__1 = 1;
00146      Treal c_b8 = 0.;
00147     
00148     /* System generated locals */
00149     integer t_dim1, t_offset, v_dim1, v_offset, i__1, i__2, i__3;
00150     Treal d__1;
00151     /* Local variables */
00152      integer i__, j;
00153      Treal vii;
00154 #define t_ref(a_1,a_2) t[(a_2)*t_dim1 + a_1]
00155 #define v_ref(a_1,a_2) v[(a_2)*v_dim1 + a_1]
00156 
00157 
00158     v_dim1 = *ldv;
00159     v_offset = 1 + v_dim1 * 1;
00160     v -= v_offset;
00161     --tau;
00162     t_dim1 = *ldt;
00163     t_offset = 1 + t_dim1 * 1;
00164     t -= t_offset;
00165 
00166     /* Function Body */
00167     if (*n == 0) {
00168         return 0;
00169     }
00170 
00171     if (template_blas_lsame(direct, "F")) {
00172         i__1 = *k;
00173         for (i__ = 1; i__ <= i__1; ++i__) {
00174             if (tau[i__] == 0.) {
00175 
00176 /*              H(i)  =  I */
00177 
00178                 i__2 = i__;
00179                 for (j = 1; j <= i__2; ++j) {
00180                     t_ref(j, i__) = 0.;
00181 /* L10: */
00182                 }
00183             } else {
00184 
00185 /*              general case */
00186 
00187                 vii = v_ref(i__, i__);
00188                 v_ref(i__, i__) = 1.;
00189                 if (template_blas_lsame(storev, "C")) {
00190 
00191 /*                 T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i) */
00192 
00193                     i__2 = *n - i__ + 1;
00194                     i__3 = i__ - 1;
00195                     d__1 = -tau[i__];
00196                     template_blas_gemv("Transpose", &i__2, &i__3, &d__1, &v_ref(i__, 1), 
00197                             ldv, &v_ref(i__, i__), &c__1, &c_b8, &t_ref(1, 
00198                             i__), &c__1);
00199                 } else {
00200 
00201 /*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)' */
00202 
00203                     i__2 = i__ - 1;
00204                     i__3 = *n - i__ + 1;
00205                     d__1 = -tau[i__];
00206                     template_blas_gemv("No transpose", &i__2, &i__3, &d__1, &v_ref(1, i__)
00207                             , ldv, &v_ref(i__, i__), ldv, &c_b8, &t_ref(1, 
00208                             i__), &c__1);
00209                 }
00210                 v_ref(i__, i__) = vii;
00211 
00212 /*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) */
00213 
00214                 i__2 = i__ - 1;
00215                 template_blas_trmv("Upper", "No transpose", "Non-unit", &i__2, &t[
00216                         t_offset], ldt, &t_ref(1, i__), &c__1);
00217                 t_ref(i__, i__) = tau[i__];
00218             }
00219 /* L20: */
00220         }
00221     } else {
00222         for (i__ = *k; i__ >= 1; --i__) {
00223             if (tau[i__] == 0.) {
00224 
00225 /*              H(i)  =  I */
00226 
00227                 i__1 = *k;
00228                 for (j = i__; j <= i__1; ++j) {
00229                     t_ref(j, i__) = 0.;
00230 /* L30: */
00231                 }
00232             } else {
00233 
00234 /*              general case */
00235 
00236                 if (i__ < *k) {
00237                     if (template_blas_lsame(storev, "C")) {
00238                         vii = v_ref(*n - *k + i__, i__);
00239                         v_ref(*n - *k + i__, i__) = 1.;
00240 
00241 /*                    T(i+1:k,i) :=   
00242                               - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i) */
00243 
00244                         i__1 = *n - *k + i__;
00245                         i__2 = *k - i__;
00246                         d__1 = -tau[i__];
00247                         template_blas_gemv("Transpose", &i__1, &i__2, &d__1, &v_ref(1, 
00248                                 i__ + 1), ldv, &v_ref(1, i__), &c__1, &c_b8, &
00249                                 t_ref(i__ + 1, i__), &c__1);
00250                         v_ref(*n - *k + i__, i__) = vii;
00251                     } else {
00252                         vii = v_ref(i__, *n - *k + i__);
00253                         v_ref(i__, *n - *k + i__) = 1.;
00254 
00255 /*                    T(i+1:k,i) :=   
00256                               - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)' */
00257 
00258                         i__1 = *k - i__;
00259                         i__2 = *n - *k + i__;
00260                         d__1 = -tau[i__];
00261                         template_blas_gemv("No transpose", &i__1, &i__2, &d__1, &v_ref(
00262                                 i__ + 1, 1), ldv, &v_ref(i__, 1), ldv, &c_b8, 
00263                                 &t_ref(i__ + 1, i__), &c__1);
00264                         v_ref(i__, *n - *k + i__) = vii;
00265                     }
00266 
00267 /*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) */
00268 
00269                     i__1 = *k - i__;
00270                     template_blas_trmv("Lower", "No transpose", "Non-unit", &i__1, &t_ref(
00271                             i__ + 1, i__ + 1), ldt, &t_ref(i__ + 1, i__), &
00272                             c__1);
00273                 }
00274                 t_ref(i__, i__) = tau[i__];
00275             }
00276 /* L40: */
00277         }
00278     }
00279     return 0;
00280 
00281 /*     End of DLARFT */
00282 
00283 } /* dlarft_ */
00284 
00285 #undef v_ref
00286 #undef t_ref
00287 
00288 
00289 #endif

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