2 #include "gromacs/utility/real.h"
4 #include "../gmx_blas.h"
5 #include "../gmx_lapack.h"
6 #include "lapack_limits.h"
10 F77_FUNC(dgesdd,DGESDD)(const char *jobz,
25 int a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
31 int idum[1], ierr, itau;
32 int minmn, wrkbl, itaup, itauq, mnthr;
35 int wntqn, wntqo, wntqs;
38 int minwrk, ldwrku, maxwrk, ldwkvt;
39 double smlnum,minval, safemin;
48 a_offset = 1 + a_dim1;
52 u_offset = 1 + u_dim1;
55 vt_offset = 1 + vt_dim1;
61 minmn = (*m < *n) ? *m : *n;
62 mnthr = (int) (minmn * 11. / 6.);
63 wntqa = (*jobz=='a' || *jobz=='A');
64 wntqs = (*jobz=='s' || *jobz=='S');
65 wntqas = wntqa || wntqs;
66 wntqn = (*jobz=='o' || *jobz=='O');
67 wntqo = (*jobz=='n' || *jobz=='N');
71 lquery = *lwork == -1;
73 if (*info == 0 && *m > 0 && *n > 0) {
79 bdspac = *n * 3 * *n + (*n << 2);
85 i__1 = wrkbl, i__2 = bdspac + *n;
86 maxwrk = (i__1 > i__2) ? i__1 : i__2;
91 i__1 = wrkbl, i__2 = *n + (*m << 5);
92 wrkbl = (i__1 > i__2) ? i__1 : i__2;
93 i__1 = wrkbl, i__2 = bdspac + *n * 3;
94 wrkbl = (i__1 > i__2) ? i__1 : i__2;
95 maxwrk = wrkbl + *n * *n;
96 minwrk = bdspac + *n * *n + *n * 3;
100 wrkbl = *n * 3 + (*m + *n*32);
102 i__1 = wrkbl, i__2 = bdspac + *n * 3;
103 maxwrk = (i__1 > i__2) ? i__1 : i__2;
104 minwrk = *n * 3 + ((*m > bdspac) ? *m : bdspac);
106 i__1 = maxwrk, i__2 = bdspac + *n * 3;
107 maxwrk = (i__1 > i__2) ? i__1 : i__2;
108 minwrk = *n * 3 + ((*m > bdspac) ? *m : bdspac);
116 bdspac = *m * 3 * *m + (*m*4);
122 i__1 = wrkbl, i__2 = bdspac + *m;
123 maxwrk = (i__1 > i__2) ? i__1 : i__2;
124 minwrk = bdspac + *m;
128 i__1 = wrkbl, i__2 = *m + (*n*32);
129 wrkbl = (i__1 > i__2) ? i__1 : i__2;
131 i__1 = wrkbl, i__2 = bdspac + *m * 3;
132 wrkbl = (i__1 > i__2) ? i__1 : i__2;
133 maxwrk = wrkbl + *m * *m;
134 minwrk = bdspac + *m * *m + *m * 3;
137 wrkbl = *m * 3 + (*m + *n*32);
139 i__1 = wrkbl, i__2 = bdspac + *m * 3;
140 maxwrk = (i__1 > i__2) ? i__1 : i__2;
141 minwrk = *m * 3 + ((*m > bdspac) ? *m : bdspac);
143 i__1 = wrkbl, i__2 = bdspac + *m * 3;
144 maxwrk = (i__1 > i__2) ? i__1 : i__2;
145 minwrk = *m * 3 + ((*m > bdspac) ? *m : bdspac);
149 work[1] = (double) maxwrk;
159 if (*m == 0 || *n == 0) {
165 eps = GMX_DOUBLE_EPS;
166 minval = GMX_DOUBLE_MIN;
167 safemin = minval / eps;
168 smlnum = sqrt(safemin) / eps;
171 bignum = 1. / smlnum;
174 anrm = F77_FUNC(dlange,DLANGE)("M", m, n, &a[a_offset], lda, dum);
176 if (anrm > 0. && anrm < smlnum) {
178 F77_FUNC(dlascl,DLASCL)("G",&c__0,&c__0,&anrm,&smlnum,m,n,&a[a_offset],lda,&ierr);
179 } else if (anrm > bignum) {
181 F77_FUNC(dlascl,DLASCL)("G",&c__0,&c__0,&anrm,&bignum,m,n,&a[a_offset],lda,&ierr);
192 i__1 = *lwork - nwork + 1;
193 F77_FUNC(dgeqrf,DGEQRF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
198 F77_FUNC(dlaset,DLASET)("L", &i__1, &i__2, &zero, &zero, &a[a_dim1 + 2],
205 i__1 = *lwork - nwork + 1;
206 F77_FUNC(dgebrd,DGEBRD)(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
207 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
210 F77_FUNC(dbdsdc,DBDSDC)("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
211 dum, idum, &work[nwork], &iwork[1], info);
217 itau = iu + ldwrku * *n;
220 i__1 = *lwork - nwork + 1;
221 F77_FUNC(dgeqrf,DGEQRF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
223 F77_FUNC(dlacpy,DLACPY)("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
225 i__1 = *lwork - nwork + 1;
226 F77_FUNC(dorgqr,DORGQR)(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork],
231 F77_FUNC(dlaset,DLASET)("L", &i__1, &i__2, &zero, &zero, &a[a_dim1 + 2],
238 i__1 = *lwork - nwork + 1;
239 F77_FUNC(dgebrd,DGEBRD)(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
240 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
242 F77_FUNC(dbdsdc,DBDSDC)("U", "I", n, &s[1], &work[ie], &work[iu], n, &vt[
243 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
246 i__1 = *lwork - nwork + 1;
247 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", n, n, n, &a[a_offset], lda, &work[
248 itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
250 i__1 = *lwork - nwork + 1;
251 F77_FUNC(dormbr,DORMBR)("P", "R", "T", n, n, n, &a[a_offset], lda, &work[
252 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
255 F77_FUNC(dgemm,DGEMM)("N", "N", m, n, n, &one, &u[u_offset], ldu, &work[iu]
256 , &ldwrku, &zero, &a[a_offset], lda);
258 F77_FUNC(dlacpy,DLACPY)("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
268 i__1 = *lwork - nwork + 1;
269 F77_FUNC(dgebrd,DGEBRD)(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
270 work[itaup], &work[nwork], &i__1, &ierr);
273 F77_FUNC(dbdsdc,DBDSDC)("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
274 dum, idum, &work[nwork], &iwork[1], info);
277 F77_FUNC(dlaset,DLASET)("F", m, m, &zero, &zero, &u[u_offset], ldu);
278 F77_FUNC(dbdsdc,DBDSDC)("U", "I", n, &s[1], &work[ie], &u[u_offset], ldu, &vt[
279 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
284 F77_FUNC(dlaset,DLASET)("F", &i__1, &i__2, &zero, &one, &u[*n + 1 + (*n +
287 i__1 = *lwork - nwork + 1;
288 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
289 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
290 i__1 = *lwork - nwork + 1;
291 F77_FUNC(dormbr,DORMBR)("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
292 itaup], &vt[vt_offset],ldvt,&work[nwork],&i__1,&ierr);
306 i__1 = *lwork - nwork + 1;
307 F77_FUNC(dgelqf,DGELQF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
312 F77_FUNC(dlaset,DLASET)("U", &i__1, &i__2, &zero, &zero, &a[(a_dim1*2) +
319 i__1 = *lwork - nwork + 1;
320 F77_FUNC(dgebrd,DGEBRD)(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
321 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
324 F77_FUNC(dbdsdc,DBDSDC)("U", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
325 dum, idum, &work[nwork], &iwork[1], info);
332 itau = ivt + ldwkvt * *m;
335 i__1 = *lwork - nwork + 1;
336 F77_FUNC(dgelqf,DGELQF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
338 F77_FUNC(dlacpy,DLACPY)("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
340 i__1 = *lwork - nwork + 1;
341 F77_FUNC(dorglq,DORGLQ)(n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[
342 nwork], &i__1, &ierr);
346 F77_FUNC(dlaset,DLASET)("U", &i__1, &i__2, &zero, &zero, &a[(a_dim1*2) +
353 i__1 = *lwork - nwork + 1;
354 F77_FUNC(dgebrd,DGEBRD)(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
355 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
357 F77_FUNC(dbdsdc,DBDSDC)("U", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &
358 work[ivt], &ldwkvt, dum, idum, &work[nwork], &iwork[1]
361 i__1 = *lwork - nwork + 1;
362 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", m, m, m, &a[a_offset], lda, &work[
363 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
364 i__1 = *lwork - nwork + 1;
365 F77_FUNC(dormbr,DORMBR)("P", "R", "T", m, m, m, &a[a_offset], lda, &work[
366 itaup], &work[ivt], &ldwkvt, &work[nwork], &i__1, &
369 F77_FUNC(dgemm,DGEMM)("N", "N", m, n, m, &one, &work[ivt], &ldwkvt, &vt[
370 vt_offset], ldvt, &zero, &a[a_offset], lda);
372 F77_FUNC(dlacpy,DLACPY)("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
383 i__1 = *lwork - nwork + 1;
384 F77_FUNC(dgebrd,DGEBRD)(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
385 work[itaup], &work[nwork], &i__1, &ierr);
388 F77_FUNC(dbdsdc,DBDSDC)("L", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
389 dum, idum, &work[nwork], &iwork[1], info);
391 F77_FUNC(dlaset,DLASET)("F", n, n, &zero, &zero, &vt[vt_offset], ldvt);
392 F77_FUNC(dbdsdc,DBDSDC)("L", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &vt[
393 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
398 F77_FUNC(dlaset,DLASET)("F", &i__1, &i__2, &zero, &one, &vt[*m + 1 + (*m +
399 1) * vt_dim1], ldvt);
401 i__1 = *lwork - nwork + 1;
402 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
403 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
404 i__1 = *lwork - nwork + 1;
405 F77_FUNC(dormbr,DORMBR)("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
406 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
416 F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
420 F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
425 work[1] = (double) maxwrk;