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(sbdsdc,SBDSDC)10 F77_FUNC(sbdsdc,SBDSDC)(const char *uplo,
11 const char *compq,
12 int *n,
13 float *d__,
14 float *e,
15 float *u,
16 int *ldu,
17 float *vt,
18 int *ldvt,
19 float *q,
20 int *iq,
21 float *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 float p, r__;
28 int z__, ic, ii, kk;
29 float cs;
30 int is, iu;
31 float sn;
32 int nm1;
33 float eps;
34 int ivt, difl, difr, ierr, perm, mlvl, sqre;
35 int poles, iuplo, nsize, start;
36 int givcol;
37 int icompq;
38 float orgnrm;
39 int givnum, givptr, qstart, smlsiz, wstart, smlszp;
40 float zero = 0.0;
41 float 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(scopy,SCOPY)(n, &d__[1], &c_1, &q[1], &c_1);
102 i__1 = *n - 1;
103 F77_FUNC(scopy,SCOPY)(&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(slartg,SLARTG)(&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(slasdq,SLASDQ)("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(slaset,SLASET)("A", n, n, &zero, &one, &u[u_offset], ldu);
131 F77_FUNC(slaset,SLASET)("A", n, n, &zero, &one, &vt[vt_offset], ldvt);
132 F77_FUNC(slasdq,SLASDQ)("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(slaset,SLASET)("A", n, n, &zero, &one, &q[iu + (qstart - 1) * *n], n);
138 F77_FUNC(slaset,SLASET)("A", n, n, &zero, &one, &q[ivt + (qstart - 1) * *n], n);
139 F77_FUNC(slasdq,SLASDQ)("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(slaset,SLASET)("A", n, n, &zero, &one, &u[u_offset], ldu);
148 F77_FUNC(slaset,SLASET)("A", n, n, &zero, &one, &vt[vt_offset], ldvt);
149 }
150
151 orgnrm = F77_FUNC(slanst,SLANST)("M", n, &d__[1], &e[1]);
152 if ( std::abs(orgnrm)<GMX_FLOAT_MIN) {
153 return;
154 }
155 F77_FUNC(slascl,SLASCL)("G", &c_0, &c_0, &orgnrm, &one, n, &c_1, &d__[1], n, &ierr);
156 F77_FUNC(slascl,SLASCL)("G", &c_0, &c_0, &orgnrm, &one, &nm1, &c_1, &e[1], &nm1, &ierr);
157
158 eps = GMX_FLOAT_EPS;
159
160 mlvl = (int) (std::log((float) (*n) / (float) (smlsiz + 1)) / std::log(2.)) + 1;
161 smlszp = smlsiz + 1;
162
163 if (icompq == 1) {
164 iu = 1;
165 ivt = smlsiz + 1;
166 difl = ivt + smlszp;
167 difr = difl + mlvl;
168 z__ = difr + (mlvl << 1);
169 ic = z__ + mlvl;
170 is = ic + 1;
171 poles = is + 1;
172 givnum = poles + (mlvl << 1);
173
174 k = 1;
175 givptr = 2;
176 perm = 3;
177 givcol = perm + mlvl;
178 }
179
180 i__1 = *n;
181 for (i__ = 1; i__ <= i__1; ++i__) {
182 if (std::abs(d__[i__]) < eps)
183 d__[i__] = (d__[i__]>0) ? eps : -eps;
184 }
185
186 start = 1;
187 sqre = 0;
188
189 i__1 = nm1;
190 for (i__ = 1; i__ <= i__1; ++i__) {
191 if (std::abs(e[i__]) < eps || i__ == nm1) {
192 if (i__ < nm1) {
193 nsize = i__ - start + 1;
194 } else if (std::abs(e[i__]) >= eps) {
195 nsize = *n - start + 1;
196 } else {
197 nsize = i__ - start + 1;
198 if (icompq == 2) {
199 u[*n + *n * u_dim1] = (d__[*n]>0) ? 1.0 : -1.0;
200 vt[*n + *n * vt_dim1] = 1.;
201 } else if (icompq == 1) {
202 q[*n + (qstart - 1) * *n] = (d__[*n]>0) ? 1.0 : -1.0;
203 q[*n + (smlsiz + qstart - 1) * *n] = 1.;
204 }
205 d__[*n] = std::abs(d__[*n]);
206 }
207 if (icompq == 2) {
208 F77_FUNC(slasd0,SLASD0)(&nsize, &sqre, &d__[start], &e[start],
209 &u[start + start * u_dim1], ldu,
210 &vt[start + start * vt_dim1],
211 ldvt, &smlsiz, &iwork[1], &work[wstart], info);
212 } else {
213 F77_FUNC(slasda,SLASDA)(&icompq, &smlsiz, &nsize, &sqre, &d__[start],
214 &e[start], &q[start + (iu + qstart - 2) * *n], n,
215 &q[start + (ivt + qstart - 2) * *n], &iq[start + k * *n],
216 &q[start + (difl + qstart - 2) * *n],
217 &q[start + (difr + qstart - 2) * *n],
218 &q[start + (z__ + qstart - 2) * *n],
219 &q[start + (poles + qstart - 2) * *n],
220 &iq[start + givptr * *n], &iq[start + givcol * *n], n,
221 &iq[start + perm * *n],
222 &q[start + (givnum + qstart - 2) * *n],
223 &q[start + (ic + qstart - 2) * *n],
224 &q[start + (is + qstart - 2) * *n], &work[wstart],
225 &iwork[1], info);
226 if (*info != 0) {
227 return;
228 }
229 }
230 start = i__ + 1;
231 }
232 }
233 F77_FUNC(slascl,SLASCL)("G", &c_0, &c_0, &one, &orgnrm, n, &c_1, &d__[1], n, &ierr);
234 L40:
235 i__1 = *n;
236 for (ii = 2; ii <= i__1; ++ii) {
237 i__ = ii - 1;
238 kk = i__;
239 p = d__[i__];
240 i__2 = *n;
241 for (j = ii; j <= i__2; ++j) {
242 if (d__[j] > p) {
243 kk = j;
244 p = d__[j];
245 }
246 }
247 if (kk != i__) {
248 d__[kk] = d__[i__];
249 d__[i__] = p;
250 if (icompq == 1) {
251 iq[i__] = kk;
252 } else if (icompq == 2) {
253 F77_FUNC(sswap,SSWAP)(n, &u[i__ * u_dim1 + 1],&c_1,&u[kk*u_dim1+1],&c_1);
254 F77_FUNC(sswap,SSWAP)(n, &vt[i__ + vt_dim1], ldvt, &vt[kk + vt_dim1], ldvt);
255 }
256 } else if (icompq == 1) {
257 iq[i__] = i__;
258 }
259 }
260 if (icompq == 1) {
261 if (iuplo == 1) {
262 iq[*n] = 1;
263 } else {
264 iq[*n] = 0;
265 }
266 }
267 if (iuplo == 2 && icompq == 2) {
268 F77_FUNC(slasr,SLASR)("L", "V", "B", n, n, &work[1], &work[*n], &u[u_offset], ldu);
269 }
270
271 return;
272 }
273