src/gromacs/linearalgebra/gmx_lapack/dormtr.c

   1 #include "../gmx_lapack.h"
   2 #include "lapack_limits.h"
   3
   4
   5 void
   6 F77_FUNC(dormtr,DORMTR)(const char *side,
   7         const char *uplo,
   8         const char *trans,
   9         int *m,
  10         int *n,
  11         double *a,
  12         int *lda,
  13         double *tau,
  14         double *c__,
  15         int *ldc,
  16         double *work,
  17         int *lwork,
  18         int *info)
  19 {
  20     int a_dim1, a_offset, c_dim1, c_offset, i__2;
  21
  22     int i1, i2, nb, mi, ni, nq, nw;
  23     int left;
  24     int iinfo;
  25     int upper;
  26     int lwkopt;
  27     int lquery;
  28
  29
  30     a_dim1 = *lda;
  31     a_offset = 1 + a_dim1;
  32     a -= a_offset;
  33     --tau;
  34     c_dim1 = *ldc;
  35     c_offset = 1 + c_dim1;
  36     c__ -= c_offset;
  37     --work;
  38
  39     *info = 0;
  40     left = (*side=='L' || *side=='l');
  41     upper = (*uplo=='U' || *uplo=='u');
  42     lquery = *lwork == -1;
  43
  44     if (left) {
  45         nq = *m;
  46         nw = *n;
  47     } else {
  48         nq = *n;
  49         nw = *m;
  50     }
  51
  52
  53     nb = DORMQL_BLOCKSIZE;
  54     lwkopt = nw * nb;
  55     work[1] = (double) lwkopt;
  56
  57     if (*info != 0) {
  58         i__2 = -(*info);
  59         return;
  60     } else if (lquery) {
  61         return;
  62     }
  63
  64     if (*m == 0 || *n == 0 || nq == 1) {
  65         work[1] = 1.;
  66         return;
  67     }
  68
  69     if (left) {
  70         mi = *m - 1;
  71         ni = *n;
  72     } else {
  73         mi = *m;
  74         ni = *n - 1;
  75     }
  76
  77     if (upper) {
  78         i__2 = nq - 1;
  79         F77_FUNC(dormql,DORMQL)(side, trans, &mi, &ni, &i__2, &a[(a_dim1 << 1) + 1], lda, &
  80                 tau[1], &c__[c_offset], ldc, &work[1], lwork, &iinfo);
  81     } else {
  82         if (left) {
  83             i1 = 2;
  84             i2 = 1;
  85         } else {
  86             i1 = 1;
  87             i2 = 2;
  88         }
  89         i__2 = nq - 1;
  90         F77_FUNC(dormqr,DORMQR)(side, trans, &mi, &ni, &i__2, &a[a_dim1 + 2], lda, &tau[1], &
  91                 c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &iinfo);
  92     }
  93     work[1] = (double) lwkopt;
  94     return;
  95
  96 }
  97
  98