template_lapack_lacon.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_LACON_HEADER
00036 #define TEMPLATE_LAPACK_LACON_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_lacon(const integer *n, Treal *v, Treal *x, 
00041         integer *isgn, Treal *est, integer *kase)
00042 {
00043 /*  -- LAPACK auxiliary 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     DLACON estimates the 1-norm of a square, real matrix A.   
00053     Reverse communication is used for evaluating matrix-vector products.   
00054 
00055     Arguments   
00056     =========   
00057 
00058     N      (input) INTEGER   
00059            The order of the matrix.  N >= 1.   
00060 
00061     V      (workspace) DOUBLE PRECISION array, dimension (N)   
00062            On the final return, V = A*W,  where  EST = norm(V)/norm(W)   
00063            (W is not returned).   
00064 
00065     X      (input/output) DOUBLE PRECISION array, dimension (N)   
00066            On an intermediate return, X should be overwritten by   
00067                  A * X,   if KASE=1,   
00068                  A' * X,  if KASE=2,   
00069            and DLACON must be re-called with all the other parameters   
00070            unchanged.   
00071 
00072     ISGN   (workspace) INTEGER array, dimension (N)   
00073 
00074     EST    (output) DOUBLE PRECISION   
00075            An estimate (a lower bound) for norm(A).   
00076 
00077     KASE   (input/output) INTEGER   
00078            On the initial call to DLACON, KASE should be 0.   
00079            On an intermediate return, KASE will be 1 or 2, indicating   
00080            whether X should be overwritten by A * X  or A' * X.   
00081            On the final return from DLACON, KASE will again be 0.   
00082 
00083     Further Details   
00084     ======= =======   
00085 
00086     Contributed by Nick Higham, University of Manchester.   
00087     Originally named SONEST, dated March 16, 1988.   
00088 
00089     Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of   
00090     a real or complex matrix, with applications to condition estimation",   
00091     ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.   
00092 
00093     =====================================================================   
00094 
00095 
00096        Parameter adjustments */
00097     /* Table of constant values */
00098      integer c__1 = 1;
00099      Treal c_b11 = 1.;
00100     
00101     /* System generated locals */
00102     integer i__1;
00103     Treal d__1;
00104     /* Builtin functions */
00105     double d_sign(Treal *, Treal *);
00106     integer i_dnnt(Treal *);
00107     /* Local variables */
00108      integer iter;
00109      Treal temp;
00110      integer jump, i__, j;
00111      integer jlast;
00112      Treal altsgn, estold;
00113 
00114 
00115     --isgn;
00116     --x;
00117     --v;
00118 
00119     /* Function Body */
00120     if (*kase == 0) {
00121         i__1 = *n;
00122         for (i__ = 1; i__ <= i__1; ++i__) {
00123             x[i__] = 1. / (Treal) (*n);
00124 /* L10: */
00125         }
00126         *kase = 1;
00127         jump = 1;
00128         return 0;
00129     }
00130 
00131     switch (jump) {
00132         case 1:  goto L20;
00133         case 2:  goto L40;
00134         case 3:  goto L70;
00135         case 4:  goto L110;
00136         case 5:  goto L140;
00137     }
00138 
00139 /*     ................ ENTRY   (JUMP = 1)   
00140        FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X. */
00141 
00142 L20:
00143     if (*n == 1) {
00144         v[1] = x[1];
00145         *est = absMACRO(v[1]);
00146 /*        ... QUIT */
00147         goto L150;
00148     }
00149     *est = template_blas_asum(n, &x[1], &c__1);
00150 
00151     i__1 = *n;
00152     for (i__ = 1; i__ <= i__1; ++i__) {
00153         x[i__] = d_sign(&c_b11, &x[i__]);
00154         isgn[i__] = i_dnnt(&x[i__]);
00155 /* L30: */
00156     }
00157     *kase = 2;
00158     jump = 2;
00159     return 0;
00160 
00161 /*     ................ ENTRY   (JUMP = 2)   
00162        FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
00163 
00164 L40:
00165     j = template_blas_idamax(n, &x[1], &c__1);
00166     iter = 2;
00167 
00168 /*     MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
00169 
00170 L50:
00171     i__1 = *n;
00172     for (i__ = 1; i__ <= i__1; ++i__) {
00173         x[i__] = 0.;
00174 /* L60: */
00175     }
00176     x[j] = 1.;
00177     *kase = 1;
00178     jump = 3;
00179     return 0;
00180 
00181 /*     ................ ENTRY   (JUMP = 3)   
00182        X HAS BEEN OVERWRITTEN BY A*X. */
00183 
00184 L70:
00185     template_blas_copy(n, &x[1], &c__1, &v[1], &c__1);
00186     estold = *est;
00187     *est = template_blas_asum(n, &v[1], &c__1);
00188     i__1 = *n;
00189     for (i__ = 1; i__ <= i__1; ++i__) {
00190         d__1 = d_sign(&c_b11, &x[i__]);
00191         if (i_dnnt(&d__1) != isgn[i__]) {
00192             goto L90;
00193         }
00194 /* L80: */
00195     }
00196 /*     REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. */
00197     goto L120;
00198 
00199 L90:
00200 /*     TEST FOR CYCLING. */
00201     if (*est <= estold) {
00202         goto L120;
00203     }
00204 
00205     i__1 = *n;
00206     for (i__ = 1; i__ <= i__1; ++i__) {
00207         x[i__] = d_sign(&c_b11, &x[i__]);
00208         isgn[i__] = i_dnnt(&x[i__]);
00209 /* L100: */
00210     }
00211     *kase = 2;
00212     jump = 4;
00213     return 0;
00214 
00215 /*     ................ ENTRY   (JUMP = 4)   
00216        X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
00217 
00218 L110:
00219     jlast = j;
00220     j = template_blas_idamax(n, &x[1], &c__1);
00221     if (x[jlast] != (d__1 = x[j], absMACRO(d__1)) && iter < 5) {
00222         ++iter;
00223         goto L50;
00224     }
00225 
00226 /*     ITERATION COMPLETE.  FINAL STAGE. */
00227 
00228 L120:
00229     altsgn = 1.;
00230     i__1 = *n;
00231     for (i__ = 1; i__ <= i__1; ++i__) {
00232         x[i__] = altsgn * ((Treal) (i__ - 1) / (Treal) (*n - 1) + 
00233                 1.);
00234         altsgn = -altsgn;
00235 /* L130: */
00236     }
00237     *kase = 1;
00238     jump = 5;
00239     return 0;
00240 
00241 /*     ................ ENTRY   (JUMP = 5)   
00242        X HAS BEEN OVERWRITTEN BY A*X. */
00243 
00244 L140:
00245     temp = template_blas_asum(n, &x[1], &c__1) / (Treal) (*n * 3) * 2.;
00246     if (temp > *est) {
00247         template_blas_copy(n, &x[1], &c__1, &v[1], &c__1);
00248         *est = temp;
00249     }
00250 
00251 L150:
00252     *kase = 0;
00253     return 0;
00254 
00255 /*     End of DLACON */
00256 
00257 } /* dlacon_ */
00258 
00259 #endif

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