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