src/gmxlib/gmx_lapack/dormbr.c

   1 #include "gmx_lapack.h"
   2 #include "lapack_limits.h"
   3
   4 void
   5 F77_FUNC(dormbr,DORMBR)(const char *vect,
   6         const char *side,
   7         const char *trans,
   8         int *m,
   9         int *n,
  10         int *k,
  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__1;
  21
  22
  23     int i1, i2, nb, mi, ni, nq, nw;
  24     int left;
  25     int iinfo;
  26     int notran;
  27     int applyq;
  28     char transt[1];
  29     int lwkopt;
  30     int lquery;
  31
  32     a_dim1 = *lda;
  33     a_offset = 1 + a_dim1;
  34     a -= a_offset;
  35     --tau;
  36     c_dim1 = *ldc;
  37     c_offset = 1 + c_dim1;
  38     c__ -= c_offset;
  39     --work;
  40     *info = 0;
  41     applyq = (*vect=='Q' || *vect=='q');
  42     left = (*side=='L' || *side=='l');
  43     notran = (*trans=='N' || *trans=='n');
  44     lquery = *lwork == -1;
  45
  46     if (left) {
  47         nq = *m;
  48         nw = *n;
  49     } else {
  50         nq = *n;
  51         nw = *m;
  52     }
  53
  54     nb = DORMQR_BLOCKSIZE;
  55     lwkopt = nw * nb;
  56     work[1] = (double) lwkopt;
  57
  58     if (*info != 0) {
  59         i__1 = -(*info);
  60         return;
  61     } else if (lquery) {
  62         return;
  63     }
  64
  65     work[1] = 1.;
  66     if (*m == 0 || *n == 0) {
  67         return;
  68     }
  69
  70     if (applyq) {
  71
  72         if (nq >= *k) {
  73
  74             F77_FUNC(dormqr,DORMQR)(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
  75                     c_offset], ldc, &work[1], lwork, &iinfo);
  76         } else if (nq > 1) {
  77
  78             if (left) {
  79                 mi = *m - 1;
  80                 ni = *n;
  81                 i1 = 2;
  82                 i2 = 1;
  83             } else {
  84                 mi = *m;
  85                 ni = *n - 1;
  86                 i1 = 1;
  87                 i2 = 2;
  88             }
  89             i__1 = nq - 1;
  90             F77_FUNC(dormqr,DORMQR)(side, trans, &mi, &ni, &i__1, &a[a_dim1 + 2], lda, &tau[1]
  91                     , &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &iinfo);
  92         }
  93     } else {
  94
  95         if (notran) {
  96             *(unsigned char *)transt = 'T';
  97         } else {
  98             *(unsigned char *)transt = 'N';
  99         }
 100         if (nq > *k) {
 101
 102             F77_FUNC(dormlq,DORMLQ)(side, transt, m, n, k, &a[a_offset], lda, &tau[1], &c__[
 103                     c_offset], ldc, &work[1], lwork, &iinfo);
 104         } else if (nq > 1) {
 105
 106             if (left) {
 107                 mi = *m - 1;
 108                 ni = *n;
 109                 i1 = 2;
 110                 i2 = 1;
 111             } else {
 112                 mi = *m;
 113                 ni = *n - 1;
 114                 i1 = 1;
 115                 i2 = 2;
 116             }
 117             i__1 = nq - 1;
 118             F77_FUNC(dormlq,DORMLQ)(side, transt, &mi, &ni, &i__1, &a[(a_dim1 << 1) + 1], lda,
 119                      &tau[1], &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &
 120                     iinfo);
 121         }
 122     }
 123     work[1] = (double) lwkopt;
 124     return;
 125
 126
 127 }
 128
 129