1 #include "../gmx_lapack.h"
2 #include "lapack_limits.h"
3
4 void
F77_FUNC(sorgbr,SORGBR)5 F77_FUNC(sorgbr,SORGBR)(const char *vect,
6 int *m,
7 int *n,
8 int *k,
9 float *a,
10 int *lda,
11 float *tau,
12 float *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(sorgqr,SORGQR)(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(sorgqr,SORGQR)(&i1,&i1,&i1,&(a[*lda+1]),lda,tau,work,lwork,&iinfo);
46 }
47 }
48 } else {
49 if(*k<*n)
50 F77_FUNC(sorglq,SORGLQ)(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(sorglq,SORGLQ)(&i1,&i1,&i1,&(a[*lda+1]),lda,tau,work,lwork,&iinfo);
63 }
64 }
65 }
66 work[0] = wrksz;
67 return;
68 }
69
70