src/gmxlib/gmx_blas/dtrsm.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 <ctype.h>
  36 #include <math.h>
  37
  38 #include <types/simple.h>
  39 #include "gmx_blas.h"
  40
  41 void
  42 F77_FUNC(dtrsm,DTRSM)(const char * side,
  43        const char * uplo,
  44        const char * transa,
  45        const char * diag,
  46        int *  m__,
  47        int *  n__,
  48        double *alpha__,
  49        double *a,
  50        int *  lda__,
  51        double *b,
  52        int *  ldb__)
  53 {
  54   const char xside  = toupper(*side);
  55   const char xuplo  = toupper(*uplo);
  56   const char xtrans = toupper(*transa);
  57   const char xdiag  = toupper(*diag);
  58   int i,j,k;
  59   double temp;
  60
  61
  62   int m = *m__;
  63   int n = *n__;
  64   int lda = *lda__;
  65   int ldb = *ldb__;
  66   double alpha = *alpha__;
  67
  68   if(n<=0)
  69     return;
  70
  71   if(fabs(alpha)<GMX_DOUBLE_MIN) {
  72     for(j=0;j<n;j++)
  73       for(i=0;i<m;i++)
  74         b[j*(ldb)+i] = 0.0;
  75     return;
  76   }
  77
  78   if(xside=='L') {
  79     /* left side */
  80     if(xtrans=='N') {
  81       /* No transpose */
  82       if(xuplo=='U') {
  83         /* upper */
  84         for(j=0;j<n;j++) {
  85           if(fabs(alpha-1.0)>GMX_DOUBLE_EPS) {
  86             for(i=0;i<m;i++)
  87               b[j*(ldb)+i] *= alpha;
  88           }
  89           for(k=m-1;k>=0;k--) {
  90             if(fabs(b[j*(ldb)+k])>GMX_DOUBLE_MIN) {
  91               if(xdiag=='N')
  92                 b[j*(ldb)+k] /= a[k*(lda)+k];
  93               for(i=0;i<k;i++)
  94                 b[j*(ldb)+i] -= b[j*(ldb)+k]*a[k*(lda)+i];
  95             }
  96           }
  97         }
  98       } else {
  99         /* lower */
 100         for(j=0;j<n;j++) {
 101           if(fabs(alpha-1.0)>GMX_DOUBLE_EPS)
 102             for(i=0;i<m;i++)
 103               b[j*(ldb)+i] *= alpha;
 104           for(k=0;k<m;k++) {
 105             if(fabs(b[j*(ldb)+k])>GMX_DOUBLE_MIN) {
 106               if(xdiag=='N')
 107                 b[j*(ldb)+k] /= a[k*(lda)+k];
 108               for(i=k+1;i<m;i++)
 109                 b[j*(ldb)+i] -= b[j*(ldb)+k]*a[k*(lda)+i];
 110             }
 111           }
 112         }
 113       }
 114     } else {
 115       /* Transpose */
 116       if(xuplo=='U') {
 117         /* upper */
 118         for(j=0;j<n;j++) {
 119           for(i=0;i<m;i++) {
 120             temp = alpha * b[j*(ldb)+i];
 121             for(k=0;k<i;k++)
 122               temp -= a[i*(lda)+k] * b[j*(ldb)+k];
 123             if(xdiag=='N')
 124                 temp /= a[i*(lda)+i];
 125             b[j*(ldb)+i] = temp;
 126           }
 127         }
 128       } else {
 129         /* lower */
 130         for(j=0;j<n;j++) {
 131           for(i=m-1;i>=0;i--) {
 132             temp = alpha * b[j*(ldb)+i];
 133             for(k=i+1;k<m;k++)
 134               temp -= a[i*(lda)+k] * b[j*(ldb)+k];
 135             if(xdiag=='N')
 136                 temp /= a[i*(lda)+i];
 137             b[j*(ldb)+i] = temp;
 138           }
 139         }
 140       }
 141     }
 142   } else {
 143     /* right side */
 144     if(xtrans=='N') {
 145       /* No transpose */
 146       if(xuplo=='U') {
 147         /* upper */
 148         for(j=0;j<n;j++) {
 149           if(fabs(alpha-1.0)>GMX_DOUBLE_EPS)
 150             for(i=0;i<m;i++)
 151               b[j*(ldb)+i] *= alpha;
 152           for(k=0;k<j;k++) {
 153             if(fabs(a[j*(lda)+k])>GMX_DOUBLE_MIN) {
 154               for(i=0;i<m;i++)
 155                 b[j*(ldb)+i] -= a[j*(lda)+k]*b[k*(ldb)+i];
 156             }
 157           }
 158           if(xdiag=='N') {
 159             temp = 1.0/a[j*(lda)+j];
 160             for(i=0;i<m;i++)
 161               b[j*(ldb)+i] *= temp;
 162           }
 163         }
 164       } else {
 165         /* lower */
 166         for(j=n-1;j>=0;j--) {
 167           if(fabs(alpha)>GMX_DOUBLE_MIN)
 168             for(i=0;i<m;i++)
 169               b[j*(ldb)+i] *= alpha;
 170           for(k=j+1;k<n;k++) {
 171             if(fabs(a[j*(lda)+k])>GMX_DOUBLE_MIN) {
 172               for(i=0;i<m;i++)
 173                 b[j*(ldb)+i] -= a[j*(lda)+k]*b[k*(ldb)+i];
 174             }
 175           }
 176           if(xdiag=='N') {
 177             temp = 1.0/a[j*(lda)+j];
 178             for(i=0;i<m;i++)
 179               b[j*(ldb)+i] *= temp;
 180           }
 181         }
 182       }
 183     } else {
 184       /* Transpose */
 185       if(xuplo=='U') {
 186         /* upper */
 187         for(k=n-1;k>=0;k--) {
 188           if(xdiag=='N') {
 189             temp = 1.0/a[k*(lda)+k];
 190             for(i=0;i<m;i++)
 191               b[k*(ldb)+i] *= temp;
 192           }
 193           for(j=0;j<k;j++) {
 194             if(fabs(a[k*(lda)+j])>GMX_DOUBLE_MIN) {
 195               temp = a[k*(lda)+j];
 196               for(i=0;i<m;i++)
 197                 b[j*(ldb)+i] -= temp * b[k*(ldb)+i];
 198             }
 199           }
 200           if(fabs(alpha-1.0)>GMX_DOUBLE_EPS)
 201             for(i=0;i<m;i++)
 202               b[k*(ldb)+i] *= alpha;
 203         }
 204       } else {
 205         /* lower */
 206         for(k=0;k<n;k++) {
 207           if(xdiag=='N') {
 208             temp = 1.0/a[k*(lda)+k];
 209             for(i=0;i<m;i++)
 210               b[k*(ldb)+i] *= temp;
 211           }
 212           for(j=k+1;j<n;j++) {
 213             if(fabs(a[k*(lda)+j])>GMX_DOUBLE_MIN) {
 214               temp = a[k*(lda)+j];
 215               for(i=0;i<m;i++)
 216                 b[j*(ldb)+i] -= temp * b[k*(ldb)+i];
 217             }
 218           }
 219           if(fabs(alpha-1.0)>GMX_DOUBLE_EPS)
 220             for(i=0;i<m;i++)
 221               b[k*(ldb)+i] *= alpha;
 222         }
 223       }
 224     }
 225   }
 226 }