00001 /* Ergo, version 3.7, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 00004 * and Anastasia Kruchinina. 00005 * 00006 * This program is free software: you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation, either version 3 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00018 * 00019 * Primary academic reference: 00020 * Ergo: An open-source program for linear-scaling electronic structure 00021 * calculations, 00022 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 00023 * Kruchinina, 00024 * SoftwareX 7, 107 (2018), 00025 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 00026 * 00027 * For further information about Ergo, see <http://www.ergoscf.org>. 00028 */ 00029 00030 /* This file belongs to the template_lapack part of the Ergo source 00031 * code. The source files in the template_lapack directory are modified 00032 * versions of files originally distributed as CLAPACK, see the 00033 * Copyright/license notice in the file template_lapack/COPYING. 00034 */ 00035 00036 00037 #ifndef TEMPLATE_LAPACK_LACON_HEADER 00038 #define TEMPLATE_LAPACK_LACON_HEADER 00039 00040 00041 template<class Treal> 00042 int template_lapack_lacon(const integer *n, Treal *v, Treal *x, 00043 integer *isgn, Treal *est, integer *kase) 00044 { 00045 /* -- LAPACK auxiliary routine (version 3.0) -- 00046 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00047 Courant Institute, Argonne National Lab, and Rice University 00048 February 29, 1992 00049 00050 00051 Purpose 00052 ======= 00053 00054 DLACON estimates the 1-norm of a square, real matrix A. 00055 Reverse communication is used for evaluating matrix-vector products. 00056 00057 Arguments 00058 ========= 00059 00060 N (input) INTEGER 00061 The order of the matrix. N >= 1. 00062 00063 V (workspace) DOUBLE PRECISION array, dimension (N) 00064 On the final return, V = A*W, where EST = norm(V)/norm(W) 00065 (W is not returned). 00066 00067 X (input/output) DOUBLE PRECISION array, dimension (N) 00068 On an intermediate return, X should be overwritten by 00069 A * X, if KASE=1, 00070 A' * X, if KASE=2, 00071 and DLACON must be re-called with all the other parameters 00072 unchanged. 00073 00074 ISGN (workspace) INTEGER array, dimension (N) 00075 00076 EST (output) DOUBLE PRECISION 00077 An estimate (a lower bound) for norm(A). 00078 00079 KASE (input/output) INTEGER 00080 On the initial call to DLACON, KASE should be 0. 00081 On an intermediate return, KASE will be 1 or 2, indicating 00082 whether X should be overwritten by A * X or A' * X. 00083 On the final return from DLACON, KASE will again be 0. 00084 00085 Further Details 00086 ======= ======= 00087 00088 Contributed by Nick Higham, University of Manchester. 00089 Originally named SONEST, dated March 16, 1988. 00090 00091 Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of 00092 a real or complex matrix, with applications to condition estimation", 00093 ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. 00094 00095 ===================================================================== 00096 00097 00098 Parameter adjustments */ 00099 /* Table of constant values */ 00100 integer c__1 = 1; 00101 Treal c_b11 = 1.; 00102 00103 /* System generated locals */ 00104 integer i__1; 00105 Treal d__1; 00106 /* Builtin functions */ 00107 double d_sign(Treal *, Treal *); 00108 integer i_dnnt(Treal *); 00109 /* Local variables */ 00110 integer iter; 00111 Treal temp; 00112 integer jump, i__, j; 00113 integer jlast; 00114 Treal altsgn, estold; 00115 00116 00117 --isgn; 00118 --x; 00119 --v; 00120 00121 /* Function Body */ 00122 if (*kase == 0) { 00123 i__1 = *n; 00124 for (i__ = 1; i__ <= i__1; ++i__) { 00125 x[i__] = 1. / (Treal) (*n); 00126 /* L10: */ 00127 } 00128 *kase = 1; 00129 jump = 1; 00130 return 0; 00131 } 00132 00133 switch (jump) { 00134 case 1: goto L20; 00135 case 2: goto L40; 00136 case 3: goto L70; 00137 case 4: goto L110; 00138 case 5: goto L140; 00139 } 00140 00141 /* ................ ENTRY (JUMP = 1) 00142 FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. */ 00143 00144 L20: 00145 if (*n == 1) { 00146 v[1] = x[1]; 00147 *est = absMACRO(v[1]); 00148 /* ... QUIT */ 00149 goto L150; 00150 } 00151 *est = template_blas_asum(n, &x[1], &c__1); 00152 00153 i__1 = *n; 00154 for (i__ = 1; i__ <= i__1; ++i__) { 00155 x[i__] = d_sign(&c_b11, &x[i__]); 00156 isgn[i__] = i_dnnt(&x[i__]); 00157 /* L30: */ 00158 } 00159 *kase = 2; 00160 jump = 2; 00161 return 0; 00162 00163 /* ................ ENTRY (JUMP = 2) 00164 FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */ 00165 00166 L40: 00167 j = template_blas_idamax(n, &x[1], &c__1); 00168 iter = 2; 00169 00170 /* MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */ 00171 00172 L50: 00173 i__1 = *n; 00174 for (i__ = 1; i__ <= i__1; ++i__) { 00175 x[i__] = 0.; 00176 /* L60: */ 00177 } 00178 x[j] = 1.; 00179 *kase = 1; 00180 jump = 3; 00181 return 0; 00182 00183 /* ................ ENTRY (JUMP = 3) 00184 X HAS BEEN OVERWRITTEN BY A*X. */ 00185 00186 L70: 00187 template_blas_copy(n, &x[1], &c__1, &v[1], &c__1); 00188 estold = *est; 00189 *est = template_blas_asum(n, &v[1], &c__1); 00190 i__1 = *n; 00191 for (i__ = 1; i__ <= i__1; ++i__) { 00192 d__1 = d_sign(&c_b11, &x[i__]); 00193 if (i_dnnt(&d__1) != isgn[i__]) { 00194 goto L90; 00195 } 00196 /* L80: */ 00197 } 00198 /* REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. */ 00199 goto L120; 00200 00201 L90: 00202 /* TEST FOR CYCLING. */ 00203 if (*est <= estold) { 00204 goto L120; 00205 } 00206 00207 i__1 = *n; 00208 for (i__ = 1; i__ <= i__1; ++i__) { 00209 x[i__] = d_sign(&c_b11, &x[i__]); 00210 isgn[i__] = i_dnnt(&x[i__]); 00211 /* L100: */ 00212 } 00213 *kase = 2; 00214 jump = 4; 00215 return 0; 00216 00217 /* ................ ENTRY (JUMP = 4) 00218 X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */ 00219 00220 L110: 00221 jlast = j; 00222 j = template_blas_idamax(n, &x[1], &c__1); 00223 if (x[jlast] != (d__1 = x[j], absMACRO(d__1)) && iter < 5) { 00224 ++iter; 00225 goto L50; 00226 } 00227 00228 /* ITERATION COMPLETE. FINAL STAGE. */ 00229 00230 L120: 00231 altsgn = 1.; 00232 i__1 = *n; 00233 for (i__ = 1; i__ <= i__1; ++i__) { 00234 x[i__] = altsgn * ((Treal) (i__ - 1) / (Treal) (*n - 1) + 00235 1.); 00236 altsgn = -altsgn; 00237 /* L130: */ 00238 } 00239 *kase = 1; 00240 jump = 5; 00241 return 0; 00242 00243 /* ................ ENTRY (JUMP = 5) 00244 X HAS BEEN OVERWRITTEN BY A*X. */ 00245 00246 L140: 00247 temp = template_blas_asum(n, &x[1], &c__1) / (Treal) (*n * 3) * 2.; 00248 if (temp > *est) { 00249 template_blas_copy(n, &x[1], &c__1, &v[1], &c__1); 00250 *est = temp; 00251 } 00252 00253 L150: 00254 *kase = 0; 00255 return 0; 00256 00257 /* End of DLACON */ 00258 00259 } /* dlacon_ */ 00260 00261 #endif