src/gmxlib/gmx_lapack/dorgbr.c

   1 #include "gmx_lapack.h"
   2 #include "lapack_limits.h"
   3
   4 void
   5 F77_FUNC(dorgbr,DORGBR)(const char *vect,
   6         int *m,
   7         int *n,
   8         int *k,
   9         double *a,
  10         int *lda,
  11         double *tau,
  12         double *work,
  13         int *lwork,
  14         int *info)
  15 {
  16   int wantq,iinfo,j,i,i1,wrksz;
  17   int mn = (*m < *n) ? *m : *n;
  18
  19   wantq = (*vect=='Q' || *vect=='q');
  20
  21   *info = 0;
  22   wrksz = mn*DORGBR_BLOCKSIZE;
  23   if(*lwork==-1) {
  24     work[0] = wrksz;
  25     return;
  26   }
  27
  28   if(*m==0 || *n==0)
  29     return;
  30
  31   if(wantq) {
  32     if(*m>=*k)
  33       F77_FUNC(dorgqr,DORGQR)(m,n,k,a,lda,tau,work,lwork,&iinfo);
  34     else {
  35       for(j=*m;j>=2;j--) {
  36         a[(j-1)*(*lda)+0] = 0.0;
  37         for(i=j+1;i<=*m;i++)
  38           a[(j-1)*(*lda)+(i-1)] = a[(j-2)*(*lda)+(i-1)];
  39       }
  40       a[0] = 1.0;
  41       for(i=2;i<=*m;i++)
  42         a[i-1] = 0.0;
  43       if(*m>1) {
  44         i1 = *m-1;
  45         F77_FUNC(dorgqr,DORGQR)(&i1,&i1,&i1,&(a[*lda+1]),lda,tau,work,lwork,&iinfo);
  46       }
  47     }
  48   } else {
  49     if(*k<*n)
  50       F77_FUNC(dorglq,DORGLQ)(m,n,k,a,lda,tau,work,lwork,&iinfo);
  51     else {
  52       a[0] = 1.0;
  53       for(i=2;i<=*m;i++)
  54         a[i-1] = 0.0;
  55       for(j=2;j<=*n;j++) {
  56         for(i=j-1;i>=2;i--)
  57           a[(j-1)*(*lda)+(i-1)] = a[(j-1)*(*lda)+(i-2)];
  58         a[(j-1)*(*lda)+0] = 0.0;
  59       }
  60       if(*n>1) {
  61         i1 = *n-1;
  62         F77_FUNC(dorglq,DORGLQ)(&i1,&i1,&i1,&(a[*lda+1]),lda,tau,work,lwork,&iinfo);
  63       }
  64     }
  65   }
  66   work[0] = wrksz;
  67   return;
  68 }
  69