template_lapack_lansy.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_LANSY_HEADER
00036 #define TEMPLATE_LAPACK_LANSY_HEADER
00037 
00038 
00039 template<class Treal>
00040 Treal template_lapack_lansy(const char *norm, const char *uplo, const integer *n, const Treal *a, const integer 
00041         *lda, Treal *work)
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        October 31, 1992   
00047 
00048 
00049     Purpose   
00050     =======   
00051 
00052     DLANSY  returns the value of the one norm,  or the Frobenius norm, or   
00053     the  infinity norm,  or the  element of  largest absolute value  of a   
00054     real symmetric matrix A.   
00055 
00056     Description   
00057     ===========   
00058 
00059     DLANSY returns the value   
00060 
00061        DLANSY = ( max(abs(A(i,j))), NORM = 'M' or 'm'   
00062                 (   
00063                 ( norm1(A),         NORM = '1', 'O' or 'o'   
00064                 (   
00065                 ( normI(A),         NORM = 'I' or 'i'   
00066                 (   
00067                 ( normF(A),         NORM = 'F', 'f', 'E' or 'e'   
00068 
00069     where  norm1  denotes the  one norm of a matrix (maximum column sum),   
00070     normI  denotes the  infinity norm  of a matrix  (maximum row sum) and   
00071     normF  denotes the  Frobenius norm of a matrix (square root of sum of   
00072     squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm.   
00073 
00074     Arguments   
00075     =========   
00076 
00077     NORM    (input) CHARACTER*1   
00078             Specifies the value to be returned in DLANSY as described   
00079             above.   
00080 
00081     UPLO    (input) CHARACTER*1   
00082             Specifies whether the upper or lower triangular part of the   
00083             symmetric matrix A is to be referenced.   
00084             = 'U':  Upper triangular part of A is referenced   
00085             = 'L':  Lower triangular part of A is referenced   
00086 
00087     N       (input) INTEGER   
00088             The order of the matrix A.  N >= 0.  When N = 0, DLANSY is   
00089             set to zero.   
00090 
00091     A       (input) DOUBLE PRECISION array, dimension (LDA,N)   
00092             The symmetric matrix A.  If UPLO = 'U', the leading n by n   
00093             upper triangular part of A contains the upper triangular part   
00094             of the matrix A, and the strictly lower triangular part of A   
00095             is not referenced.  If UPLO = 'L', the leading n by n lower   
00096             triangular part of A contains the lower triangular part of   
00097             the matrix A, and the strictly upper triangular part of A is   
00098             not referenced.   
00099 
00100     LDA     (input) INTEGER   
00101             The leading dimension of the array A.  LDA >= max(N,1).   
00102 
00103     WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK),   
00104             where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,   
00105             WORK is not referenced.   
00106 
00107    =====================================================================   
00108 
00109 
00110        Parameter adjustments */
00111     /* Table of constant values */
00112      integer c__1 = 1;
00113     
00114     /* System generated locals */
00115     integer a_dim1, a_offset, i__1, i__2;
00116     Treal ret_val, d__1, d__2, d__3;
00117     /* Local variables */
00118      Treal absa;
00119      integer i__, j;
00120      Treal scale;
00121      Treal value;
00122      Treal sum;
00123 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00124 
00125 
00126     a_dim1 = *lda;
00127     a_offset = 1 + a_dim1 * 1;
00128     a -= a_offset;
00129     --work;
00130 
00131     /* Initialization added by Elias to get rid of compiler warnings. */
00132     value = 0;
00133     /* Function Body */
00134     if (*n == 0) {
00135         value = 0.;
00136     } else if (template_blas_lsame(norm, "M")) {
00137 
00138 /*        Find max(abs(A(i,j))). */
00139 
00140         value = 0.;
00141         if (template_blas_lsame(uplo, "U")) {
00142             i__1 = *n;
00143             for (j = 1; j <= i__1; ++j) {
00144                 i__2 = j;
00145                 for (i__ = 1; i__ <= i__2; ++i__) {
00146 /* Computing MAX */
00147                     d__2 = value, d__3 = (d__1 = a_ref(i__, j), absMACRO(d__1));
00148                     value = maxMACRO(d__2,d__3);
00149 /* L10: */
00150                 }
00151 /* L20: */
00152             }
00153         } else {
00154             i__1 = *n;
00155             for (j = 1; j <= i__1; ++j) {
00156                 i__2 = *n;
00157                 for (i__ = j; i__ <= i__2; ++i__) {
00158 /* Computing MAX */
00159                     d__2 = value, d__3 = (d__1 = a_ref(i__, j), absMACRO(d__1));
00160                     value = maxMACRO(d__2,d__3);
00161 /* L30: */
00162                 }
00163 /* L40: */
00164             }
00165         }
00166     } else if (template_blas_lsame(norm, "I") || template_blas_lsame(norm, "O") || *(unsigned char *)norm == '1') {
00167 
00168 /*        Find normI(A) ( = norm1(A), since A is symmetric). */
00169 
00170         value = 0.;
00171         if (template_blas_lsame(uplo, "U")) {
00172             i__1 = *n;
00173             for (j = 1; j <= i__1; ++j) {
00174                 sum = 0.;
00175                 i__2 = j - 1;
00176                 for (i__ = 1; i__ <= i__2; ++i__) {
00177                     absa = (d__1 = a_ref(i__, j), absMACRO(d__1));
00178                     sum += absa;
00179                     work[i__] += absa;
00180 /* L50: */
00181                 }
00182                 work[j] = sum + (d__1 = a_ref(j, j), absMACRO(d__1));
00183 /* L60: */
00184             }
00185             i__1 = *n;
00186             for (i__ = 1; i__ <= i__1; ++i__) {
00187 /* Computing MAX */
00188                 d__1 = value, d__2 = work[i__];
00189                 value = maxMACRO(d__1,d__2);
00190 /* L70: */
00191             }
00192         } else {
00193             i__1 = *n;
00194             for (i__ = 1; i__ <= i__1; ++i__) {
00195                 work[i__] = 0.;
00196 /* L80: */
00197             }
00198             i__1 = *n;
00199             for (j = 1; j <= i__1; ++j) {
00200                 sum = work[j] + (d__1 = a_ref(j, j), absMACRO(d__1));
00201                 i__2 = *n;
00202                 for (i__ = j + 1; i__ <= i__2; ++i__) {
00203                     absa = (d__1 = a_ref(i__, j), absMACRO(d__1));
00204                     sum += absa;
00205                     work[i__] += absa;
00206 /* L90: */
00207                 }
00208                 value = maxMACRO(value,sum);
00209 /* L100: */
00210             }
00211         }
00212     } else if (template_blas_lsame(norm, "F") || template_blas_lsame(norm, "E")) {
00213 
00214 /*        Find normF(A). */
00215 
00216         scale = 0.;
00217         sum = 1.;
00218         if (template_blas_lsame(uplo, "U")) {
00219             i__1 = *n;
00220             for (j = 2; j <= i__1; ++j) {
00221                 i__2 = j - 1;
00222                 template_lapack_lassq(&i__2, &a_ref(1, j), &c__1, &scale, &sum);
00223 /* L110: */
00224             }
00225         } else {
00226             i__1 = *n - 1;
00227             for (j = 1; j <= i__1; ++j) {
00228                 i__2 = *n - j;
00229                 template_lapack_lassq(&i__2, &a_ref(j + 1, j), &c__1, &scale, &sum);
00230 /* L120: */
00231             }
00232         }
00233         sum *= 2;
00234         i__1 = *lda + 1;
00235         template_lapack_lassq(n, &a[a_offset], &i__1, &scale, &sum);
00236         value = scale * template_blas_sqrt(sum);
00237     }
00238 
00239     ret_val = value;
00240     return ret_val;
00241 
00242 /*     End of DLANSY */
00243 
00244 } /* dlansy_ */
00245 
00246 #undef a_ref
00247 
00248 
00249 #endif

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