src/gmxlib/gmx_lapack/dtrtri.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, 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 #include "gmx_blas.h"
  37 #include "gmx_lapack.h"
  38 #include "lapack_limits.h"
  39
  40 #include <types/simple.h>
  41
  42 void
  43 F77_FUNC(dtrtri,DTRTRI)(const char *uplo,
  44         const char *diag,
  45         int *n,
  46         double *a,
  47         int *lda,
  48         int *info)
  49 {
  50     int a_dim1, a_offset, i__1, i__3, i__4, i__5;
  51     int j, jb, nb, nn;
  52     double c_b18 = 1.;
  53     double c_b22 = -1.;
  54
  55     int upper;
  56     int nounit;
  57
  58     a_dim1 = *lda;
  59     a_offset = 1 + a_dim1;
  60     a -= a_offset;
  61
  62     *info = 0;
  63     upper = (*uplo=='U' || *uplo=='u');
  64     nounit = (*diag=='N' || *diag=='n');
  65
  66     if (*info != 0) {
  67         i__1 = -(*info);
  68         return;
  69     }
  70
  71     if (*n == 0) {
  72         return;
  73     }
  74
  75     if (nounit) {
  76         i__1 = *n;
  77         for (*info = 1; *info <= i__1; ++(*info)) {
  78             if (fabs(a[*info + *info * a_dim1])<GMX_DOUBLE_MIN) {
  79                 return;
  80             }
  81         }
  82         *info = 0;
  83     }
  84
  85     nb = DTRTRI_BLOCKSIZE;
  86     if (nb <= 1 || nb >= *n) {
  87
  88         F77_FUNC(dtrti2,DTRTI2)(uplo, diag, n, &a[a_offset], lda, info);
  89     } else {
  90
  91         if (upper) {
  92
  93             i__1 = *n;
  94             i__3 = nb;
  95             for (j = 1; i__3 < 0 ? j >= i__1 : j <= i__1; j += i__3) {
  96                 i__4 = nb, i__5 = *n - j + 1;
  97                 jb = (i__4<i__5) ? i__4 : i__5;
  98
  99                 i__4 = j - 1;
 100                 F77_FUNC(dtrmm,DTRMM)("Left", "Upper", "No transpose", diag, &i__4, &jb, &
 101                         c_b18, &a[a_offset], lda, &a[j * a_dim1 + 1], lda);
 102                 i__4 = j - 1;
 103                 F77_FUNC(dtrsm,DTRSM)("Right", "Upper", "No transpose", diag, &i__4, &jb, &
 104                         c_b22, &a[j + j * a_dim1], lda, &a[j * a_dim1 + 1],
 105                         lda);
 106
 107                 F77_FUNC(dtrti2,DTRTI2)("Upper", diag, &jb, &a[j + j * a_dim1], lda, info);
 108             }
 109         } else {
 110
 111             nn = (*n - 1) / nb * nb + 1;
 112             i__3 = -nb;
 113             for (j = nn; i__3 < 0 ? j >= 1 : j <= 1; j += i__3) {
 114                 i__1 = nb, i__4 = *n - j + 1;
 115                 jb = (i__1<i__4) ? i__1 : i__4;
 116                 if (j + jb <= *n) {
 117
 118                     i__1 = *n - j - jb + 1;
 119                     F77_FUNC(dtrmm,DTRMM)("Left", "Lower", "No transpose", diag, &i__1, &jb,
 120                             &c_b18, &a[j + jb + (j + jb) * a_dim1], lda, &a[j
 121                             + jb + j * a_dim1], lda);
 122                     i__1 = *n - j - jb + 1;
 123                     F77_FUNC(dtrsm,DTRSM)("Right", "Lower", "No transpose", diag, &i__1, &jb,
 124                              &c_b22, &a[j + j * a_dim1], lda, &a[j + jb + j *
 125                             a_dim1], lda);
 126                 }
 127
 128                 F77_FUNC(dtrti2,DTRTI2)("Lower", diag, &jb, &a[j + j * a_dim1], lda, info);
 129             }
 130         }
 131     }
 132     return;
 133 }
 134
 135