src/gmxlib/gmx_lapack/dlansy.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include <math.h>
  36
  37
  38 #include "gmx_lapack.h"
  39
  40 double
  41 F77_FUNC(dlansy,DLANSY)(const char *norm, const char *uplo, int *n, double *a, int
  42         *lda, double *work)
  43 {
  44     /* System generated locals */
  45     int a_dim1, a_offset, i__1, i__2;
  46     double ret_val, d__1, d__2, d__3;
  47     int c__1 = 1;
  48
  49     /* Local variables */
  50     int i__, j;
  51     double sum, absa, scale;
  52     double value =0.0;
  53
  54     a_dim1 = *lda;
  55     a_offset = 1 + a_dim1;
  56     a -= a_offset;
  57     --work;
  58
  59     if (*n == 0) {
  60         value = 0.;
  61     } else if (*norm=='M' || *norm=='m') {
  62
  63         value = 0.;
  64         if (*uplo=='U' || *uplo=='u') {
  65             i__1 = *n;
  66             for (j = 1; j <= i__1; ++j) {
  67                 i__2 = j;
  68                 for (i__ = 1; i__ <= i__2; ++i__) {
  69                   d__2 = value;
  70                   d__3 = fabs(a[i__ + j * a_dim1]);
  71                   value = (d__2>d__3) ? d__2 : d__3;
  72                 }
  73             }
  74         } else {
  75             i__1 = *n;
  76             for (j = 1; j <= i__1; ++j) {
  77                 i__2 = *n;
  78                 for (i__ = j; i__ <= i__2; ++i__) {
  79                   d__2 = value;
  80                   d__3 = fabs(a[i__ + j * a_dim1]);
  81                     value =  (d__2>d__3) ? d__2 : d__3;
  82                 }
  83             }
  84         }
  85     } else if (*norm=='I' || *norm=='i' || *norm=='O' || *norm=='o' || *norm=='1') {
  86
  87         value = 0.;
  88         if (*uplo=='U' || *uplo=='u') {
  89             i__1 = *n;
  90             for (j = 1; j <= i__1; ++j) {
  91                 sum = 0.;
  92                 i__2 = j - 1;
  93                 for (i__ = 1; i__ <= i__2; ++i__) {
  94                     absa = fabs(a[i__ + j * a_dim1]);
  95                     sum += absa;
  96                     work[i__] += absa;
  97                 }
  98                 work[j] = sum + fabs(a[j + j * a_dim1]);
  99             }
 100             i__1 = *n;
 101             for (i__ = 1; i__ <= i__1; ++i__) {
 102                 d__1 = value, d__2 = work[i__];
 103                 value =  (d__1>d__2) ? d__1 : d__2;
 104             }
 105         } else {
 106             i__1 = *n;
 107             for (i__ = 1; i__ <= i__1; ++i__) {
 108                 work[i__] = 0.;
 109             }
 110             i__1 = *n;
 111             for (j = 1; j <= i__1; ++j) {
 112                 sum = work[j] + fabs(a[j + j * a_dim1]);
 113                 i__2 = *n;
 114                 for (i__ = j + 1; i__ <= i__2; ++i__) {
 115                     absa = fabs(a[i__ + j * a_dim1]);
 116                     sum += absa;
 117                     work[i__] += absa;
 118                 }
 119                 if(sum>value)
 120                   value = sum;
 121             }
 122         }
 123     } else if (*norm=='F' || *norm=='f' || *norm=='E' || *norm=='e') {
 124
 125         scale = 0.;
 126         sum = 1.;
 127         if (*uplo=='U' || *uplo=='u') {
 128             i__1 = *n;
 129             for (j = 2; j <= i__1; ++j) {
 130                 i__2 = j - 1;
 131                 F77_FUNC(dlassq,DLASSQ)(&i__2, &a[j * a_dim1 + 1], &c__1, &scale, &sum);
 132             }
 133         } else {
 134             i__1 = *n - 1;
 135             for (j = 1; j <= i__1; ++j) {
 136                 i__2 = *n - j;
 137                 F77_FUNC(dlassq,DLASSQ)(&i__2, &a[j + 1 + j * a_dim1], &c__1, &scale, &sum);
 138             }
 139         }
 140         sum *= 2;
 141         i__1 = *lda + 1;
 142         F77_FUNC(dlassq,DLASSQ)(n, &a[a_offset], &i__1, &scale, &sum);
 143         value = scale * sqrt(sum);
 144     }
 145
 146     ret_val = value;
 147     return ret_val;
 148 }
 149
 150