1 #include <cctype>
2 #include <cmath>
3 #include "../gmx_blas.h"
4 #include "../gmx_lapack.h"
5 #include "lapack_limits.h"
6 
7 #include "gromacs/utility/real.h"
8 
9 void
F77_FUNC(dbdsdc,DBDSDC)10 F77_FUNC(dbdsdc,DBDSDC)(const char *uplo,
11 	const char *compq,
12 	int *n,
13 	double *d__,
14 	double *e,
15 	double *u,
16 	int *ldu,
17 	double *vt,
18 	int *ldvt,
19 	double *q,
20 	int *iq,
21 	double *work,
22 	int *iwork,
23 	int *info)
24 {
25     int u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
26     int i__, j, k;
27     double p, r__;
28     int z__, ic, ii, kk;
29     double cs;
30     int is, iu;
31     double sn;
32     int nm1;
33     double eps;
34     int ivt, difl, difr, ierr, perm, mlvl, sqre;
35     int poles, iuplo, nsize, start;
36     int givcol;
37     int icompq;
38     double orgnrm;
39     int givnum, givptr, qstart, smlsiz, wstart, smlszp;
40     double zero = 0.0;
41     double one = 1.0;
42     int c_0 = 0;
43     int c_1 = 1;
44 
45     --d__;
46     --e;
47     u_dim1 = *ldu;
48     u_offset = 1 + u_dim1;
49     u -= u_offset;
50     vt_dim1 = *ldvt;
51     vt_offset = 1 + vt_dim1;
52     vt -= vt_offset;
53     --q;
54     --iq;
55     --work;
56     --iwork;
57 
58     k = iu = z__ = ic = is = ivt = difl = difr = perm = 0;
59     poles = givnum = givptr = givcol = 0;
60 
61     smlsiz = DBDSDC_SMALLSIZE;
62     *info = 0;
63 
64     iuplo = (*uplo=='U' || *uplo=='u') ? 1 : 2;
65 
66     switch(*compq) {
67     case 'n':
68     case 'N':
69       icompq = 0;
70       break;
71     case 'p':
72     case 'P':
73       icompq = 1;
74       break;
75     case 'i':
76     case 'I':
77       icompq = 2;
78       break;
79     default:
80       return;
81     }
82 
83     if (*n <= 0)
84 	return;
85 
86     if (*n == 1) {
87 	if (icompq == 1) {
88 	  q[1] = (d__[1]>0) ? 1.0 : -1.0;
89 	  q[smlsiz * *n + 1] = 1.;
90 	} else if (icompq == 2) {
91 	  u[u_dim1 + 1] = (d__[1]>0) ? 1.0 : -1.0;
92 	  vt[vt_dim1 + 1] = 1.;
93 	}
94 	d__[1] = std::abs(d__[1]);
95 	return;
96     }
97     nm1 = *n - 1;
98     wstart = 1;
99     qstart = 3;
100     if (icompq == 1) {
101 	F77_FUNC(dcopy,DCOPY)(n, &d__[1], &c_1, &q[1], &c_1);
102 	i__1 = *n - 1;
103 	F77_FUNC(dcopy,DCOPY)(&i__1, &e[1], &c_1, &q[*n + 1], &c_1);
104     }
105     if (iuplo == 2) {
106 	qstart = 5;
107 	wstart = (*n << 1) - 1;
108 	i__1 = *n - 1;
109 	for (i__ = 1; i__ <= i__1; ++i__) {
110 	    F77_FUNC(dlartg,DLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
111 	    d__[i__] = r__;
112 	    e[i__] = sn * d__[i__ + 1];
113 	    d__[i__ + 1] = cs * d__[i__ + 1];
114 	    if (icompq == 1) {
115 		q[i__ + (*n << 1)] = cs;
116 		q[i__ + *n * 3] = sn;
117 	    } else if (icompq == 2) {
118 		work[i__] = cs;
119 		work[nm1 + i__] = -sn;
120 	    }
121 	}
122     }
123     if (icompq == 0) {
124       F77_FUNC(dlasdq,DLASDQ)("U",&c_0,n,&c_0,&c_0,&c_0,&d__[1],&e[1],&vt[vt_offset],ldvt,
125 	      &u[u_offset], ldu, &u[u_offset], ldu, &work[wstart], info);
126 	goto L40;
127     }
128     if (*n <= smlsiz) {
129 	if (icompq == 2) {
130 	    F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &u[u_offset], ldu);
131 	    F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &vt[vt_offset], ldvt);
132 	    F77_FUNC(dlasdq,DLASDQ)("U",&c_0,n,n,n,&c_0,&d__[1],&e[1],&vt[vt_offset],ldvt,
133 		    &u[u_offset],ldu,&u[u_offset],ldu,&work[wstart],info);
134 	} else if (icompq == 1) {
135 	    iu = 1;
136 	    ivt = iu + *n;
137 	    F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &q[iu + (qstart - 1) * *n], n);
138 	    F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &q[ivt + (qstart - 1) * *n], n);
139 	    F77_FUNC(dlasdq,DLASDQ)("U", &c_0, n, n, n, &c_0, &d__[1], &e[1],
140 		    &q[ivt + (qstart - 1) * *n], n, &q[iu + (qstart - 1) * *n],
141 		    n, &q[iu + (qstart - 1) * *n], n, &work[wstart], info);
142 	}
143 	goto L40;
144     }
145 
146     if (icompq == 2) {
147 	F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &u[u_offset], ldu);
148 	F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &vt[vt_offset], ldvt);
149     }
150 
151     orgnrm = F77_FUNC(dlanst,DLANST)("M", n, &d__[1], &e[1]);
152     if ( std::abs(orgnrm)<GMX_DOUBLE_MIN) {
153 	return;
154     }
155     F77_FUNC(dlascl,DLASCL)("G", &c_0, &c_0, &orgnrm, &one, n, &c_1, &d__[1], n, &ierr);
156     F77_FUNC(dlascl,DLASCL)("G", &c_0, &c_0, &orgnrm, &one, &nm1, &c_1, &e[1], &nm1, &ierr);
157 
158     eps = GMX_DOUBLE_EPS;
159 
160     mlvl = (int) (std::log((double) (*n) / (double) (smlsiz + 1)) /
161 		  std::log(2.)) + 1;
162     smlszp = smlsiz + 1;
163 
164     if (icompq == 1) {
165 	iu = 1;
166 	ivt = smlsiz + 1;
167 	difl = ivt + smlszp;
168 	difr = difl + mlvl;
169 	z__ = difr + (mlvl << 1);
170 	ic = z__ + mlvl;
171 	is = ic + 1;
172 	poles = is + 1;
173 	givnum = poles + (mlvl << 1);
174 
175 	k = 1;
176 	givptr = 2;
177 	perm = 3;
178 	givcol = perm + mlvl;
179     }
180 
181     i__1 = *n;
182     for (i__ = 1; i__ <= i__1; ++i__) {
183 	if (std::abs(d__[i__]) < eps)
184 	    d__[i__] = (d__[i__]>0) ? eps : -eps;
185     }
186 
187     start = 1;
188     sqre = 0;
189 
190     i__1 = nm1;
191     for (i__ = 1; i__ <= i__1; ++i__) {
192 	if (std::abs(e[i__]) < eps || i__ == nm1) {
193 	    if (i__ < nm1) {
194 		nsize = i__ - start + 1;
195 	    } else if (std::abs(e[i__]) >= eps) {
196 		nsize = *n - start + 1;
197 	    } else {
198 		nsize = i__ - start + 1;
199 		if (icompq == 2) {
200 		    u[*n + *n * u_dim1] = (d__[*n]>0) ? 1.0 : -1.0;
201 		    vt[*n + *n * vt_dim1] = 1.;
202 		} else if (icompq == 1) {
203 		    q[*n + (qstart - 1) * *n] = (d__[*n]>0) ? 1.0 : -1.0;
204 		    q[*n + (smlsiz + qstart - 1) * *n] = 1.;
205 		}
206 		d__[*n] = std::abs(d__[*n]);
207 	    }
208 	    if (icompq == 2) {
209 		F77_FUNC(dlasd0,DLASD0)(&nsize, &sqre, &d__[start], &e[start],
210 			&u[start + start * u_dim1], ldu,
211 			&vt[start + start * vt_dim1],
212 			ldvt, &smlsiz, &iwork[1], &work[wstart], info);
213 	    } else {
214 		F77_FUNC(dlasda,DLASDA)(&icompq, &smlsiz, &nsize, &sqre, &d__[start],
215 			&e[start], &q[start + (iu + qstart - 2) * *n], n,
216 			&q[start + (ivt + qstart - 2) * *n], &iq[start + k * *n],
217 			&q[start + (difl + qstart - 2) * *n],
218 			&q[start + (difr + qstart - 2) * *n],
219 			&q[start + (z__ + qstart - 2) * *n],
220 			&q[start + (poles + qstart - 2) * *n],
221 			&iq[start + givptr * *n], &iq[start + givcol * *n], n,
222 			&iq[start + perm * *n],
223 			&q[start + (givnum + qstart - 2) * *n],
224 			&q[start + (ic + qstart - 2) * *n],
225 			&q[start + (is + qstart - 2) * *n], &work[wstart],
226 			&iwork[1], info);
227 		if (*info != 0) {
228 		    return;
229 		}
230 	    }
231 	    start = i__ + 1;
232 	}
233     }
234     F77_FUNC(dlascl,DLASCL)("G", &c_0, &c_0, &one, &orgnrm, n, &c_1, &d__[1], n, &ierr);
235 L40:
236     i__1 = *n;
237     for (ii = 2; ii <= i__1; ++ii) {
238 	i__ = ii - 1;
239 	kk = i__;
240 	p = d__[i__];
241 	i__2 = *n;
242 	for (j = ii; j <= i__2; ++j) {
243 	    if (d__[j] > p) {
244 		kk = j;
245 		p = d__[j];
246 	    }
247 	}
248 	if (kk != i__) {
249 	    d__[kk] = d__[i__];
250 	    d__[i__] = p;
251 	    if (icompq == 1) {
252 		iq[i__] = kk;
253 	    } else if (icompq == 2) {
254 		F77_FUNC(dswap,DSWAP)(n, &u[i__ * u_dim1 + 1],&c_1,&u[kk*u_dim1+1],&c_1);
255 		F77_FUNC(dswap,DSWAP)(n, &vt[i__ + vt_dim1], ldvt, &vt[kk + vt_dim1], ldvt);
256 	    }
257 	} else if (icompq == 1) {
258 	    iq[i__] = i__;
259 	}
260     }
261     if (icompq == 1) {
262 	if (iuplo == 1) {
263 	    iq[*n] = 1;
264 	} else {
265 	    iq[*n] = 0;
266 	}
267     }
268     if (iuplo == 2 && icompq == 2) {
269 	F77_FUNC(dlasr,DLASR)("L", "V", "B", n, n, &work[1], &work[*n], &u[u_offset], ldu);
270     }
271 
272     return;
273 }
274