C_LAPACK: Fixes to make it compile with MSVC (#3605)
[platform/upstream/openblas.git] / lapack-netlib / SRC / cuncsd.c
1 #include <math.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <stdio.h>
5 #include <complex.h>
6 #ifdef complex
7 #undef complex
8 #endif
9 #ifdef I
10 #undef I
11 #endif
12
13 #if defined(_WIN64)
14 typedef long long BLASLONG;
15 typedef unsigned long long BLASULONG;
16 #else
17 typedef long BLASLONG;
18 typedef unsigned long BLASULONG;
19 #endif
20
21 #ifdef LAPACK_ILP64
22 typedef BLASLONG blasint;
23 #if defined(_WIN64)
24 #define blasabs(x) llabs(x)
25 #else
26 #define blasabs(x) labs(x)
27 #endif
28 #else
29 typedef int blasint;
30 #define blasabs(x) abs(x)
31 #endif
32
33 typedef blasint integer;
34
35 typedef unsigned int uinteger;
36 typedef char *address;
37 typedef short int shortint;
38 typedef float real;
39 typedef double doublereal;
40 typedef struct { real r, i; } complex;
41 typedef struct { doublereal r, i; } doublecomplex;
42 #ifdef _MSC_VER
43 static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
44 static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
45 static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
46 static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
47 #else
48 static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
49 static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
50 static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
51 static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
52 #endif
53 #define pCf(z) (*_pCf(z))
54 #define pCd(z) (*_pCd(z))
55 typedef int logical;
56 typedef short int shortlogical;
57 typedef char logical1;
58 typedef char integer1;
59
60 #define TRUE_ (1)
61 #define FALSE_ (0)
62
63 /* Extern is for use with -E */
64 #ifndef Extern
65 #define Extern extern
66 #endif
67
68 /* I/O stuff */
69
70 typedef int flag;
71 typedef int ftnlen;
72 typedef int ftnint;
73
74 /*external read, write*/
75 typedef struct
76 {       flag cierr;
77         ftnint ciunit;
78         flag ciend;
79         char *cifmt;
80         ftnint cirec;
81 } cilist;
82
83 /*internal read, write*/
84 typedef struct
85 {       flag icierr;
86         char *iciunit;
87         flag iciend;
88         char *icifmt;
89         ftnint icirlen;
90         ftnint icirnum;
91 } icilist;
92
93 /*open*/
94 typedef struct
95 {       flag oerr;
96         ftnint ounit;
97         char *ofnm;
98         ftnlen ofnmlen;
99         char *osta;
100         char *oacc;
101         char *ofm;
102         ftnint orl;
103         char *oblnk;
104 } olist;
105
106 /*close*/
107 typedef struct
108 {       flag cerr;
109         ftnint cunit;
110         char *csta;
111 } cllist;
112
113 /*rewind, backspace, endfile*/
114 typedef struct
115 {       flag aerr;
116         ftnint aunit;
117 } alist;
118
119 /* inquire */
120 typedef struct
121 {       flag inerr;
122         ftnint inunit;
123         char *infile;
124         ftnlen infilen;
125         ftnint  *inex;  /*parameters in standard's order*/
126         ftnint  *inopen;
127         ftnint  *innum;
128         ftnint  *innamed;
129         char    *inname;
130         ftnlen  innamlen;
131         char    *inacc;
132         ftnlen  inacclen;
133         char    *inseq;
134         ftnlen  inseqlen;
135         char    *indir;
136         ftnlen  indirlen;
137         char    *infmt;
138         ftnlen  infmtlen;
139         char    *inform;
140         ftnint  informlen;
141         char    *inunf;
142         ftnlen  inunflen;
143         ftnint  *inrecl;
144         ftnint  *innrec;
145         char    *inblank;
146         ftnlen  inblanklen;
147 } inlist;
148
149 #define VOID void
150
151 union Multitype {       /* for multiple entry points */
152         integer1 g;
153         shortint h;
154         integer i;
155         /* longint j; */
156         real r;
157         doublereal d;
158         complex c;
159         doublecomplex z;
160         };
161
162 typedef union Multitype Multitype;
163
164 struct Vardesc {        /* for Namelist */
165         char *name;
166         char *addr;
167         ftnlen *dims;
168         int  type;
169         };
170 typedef struct Vardesc Vardesc;
171
172 struct Namelist {
173         char *name;
174         Vardesc **vars;
175         int nvars;
176         };
177 typedef struct Namelist Namelist;
178
179 #define abs(x) ((x) >= 0 ? (x) : -(x))
180 #define dabs(x) (fabs(x))
181 #define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
182 #define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
183 #define dmin(a,b) (f2cmin(a,b))
184 #define dmax(a,b) (f2cmax(a,b))
185 #define bit_test(a,b)   ((a) >> (b) & 1)
186 #define bit_clear(a,b)  ((a) & ~((uinteger)1 << (b)))
187 #define bit_set(a,b)    ((a) |  ((uinteger)1 << (b)))
188
189 #define abort_() { sig_die("Fortran abort routine called", 1); }
190 #define c_abs(z) (cabsf(Cf(z)))
191 #define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
192 #ifdef _MSC_VER
193 #define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
194 #define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/df(b)._Val[1]);}
195 #else
196 #define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
197 #define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
198 #endif
199 #define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
200 #define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
201 #define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
202 //#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
203 #define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
204 #define d_abs(x) (fabs(*(x)))
205 #define d_acos(x) (acos(*(x)))
206 #define d_asin(x) (asin(*(x)))
207 #define d_atan(x) (atan(*(x)))
208 #define d_atn2(x, y) (atan2(*(x),*(y)))
209 #define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
210 #define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
211 #define d_cos(x) (cos(*(x)))
212 #define d_cosh(x) (cosh(*(x)))
213 #define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
214 #define d_exp(x) (exp(*(x)))
215 #define d_imag(z) (cimag(Cd(z)))
216 #define r_imag(z) (cimagf(Cf(z)))
217 #define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
218 #define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
219 #define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
220 #define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
221 #define d_log(x) (log(*(x)))
222 #define d_mod(x, y) (fmod(*(x), *(y)))
223 #define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
224 #define d_nint(x) u_nint(*(x))
225 #define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
226 #define d_sign(a,b) u_sign(*(a),*(b))
227 #define r_sign(a,b) u_sign(*(a),*(b))
228 #define d_sin(x) (sin(*(x)))
229 #define d_sinh(x) (sinh(*(x)))
230 #define d_sqrt(x) (sqrt(*(x)))
231 #define d_tan(x) (tan(*(x)))
232 #define d_tanh(x) (tanh(*(x)))
233 #define i_abs(x) abs(*(x))
234 #define i_dnnt(x) ((integer)u_nint(*(x)))
235 #define i_len(s, n) (n)
236 #define i_nint(x) ((integer)u_nint(*(x)))
237 #define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
238 #define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
239 #define pow_si(B,E) spow_ui(*(B),*(E))
240 #define pow_ri(B,E) spow_ui(*(B),*(E))
241 #define pow_di(B,E) dpow_ui(*(B),*(E))
242 #define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
243 #define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
244 #define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
245 #define s_cat(lpp, rpp, rnp, np, llp) {         ftnlen i, nc, ll; char *f__rp, *lp;     ll = (llp); lp = (lpp);         for(i=0; i < (int)*(np); ++i) {                 nc = ll;                if((rnp)[i] < nc) nc = (rnp)[i];                ll -= nc;               f__rp = (rpp)[i];               while(--nc >= 0) *lp++ = *(f__rp)++;         }  while(--ll >= 0) *lp++ = ' '; }
246 #define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
247 #define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
248 #define sig_die(s, kill) { exit(1); }
249 #define s_stop(s, n) {exit(0);}
250 static char junk[] = "\n@(#)LIBF77 VERSION 19990503\n";
251 #define z_abs(z) (cabs(Cd(z)))
252 #define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
253 #define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
254 #define myexit_() break;
255 #define mycycle() continue;
256 #define myceiling(w) {ceil(w)}
257 #define myhuge(w) {HUGE_VAL}
258 //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
259 #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)}
260
261 /* procedure parameter types for -A and -C++ */
262
263 #define F2C_proc_par_types 1
264 #ifdef __cplusplus
265 typedef logical (*L_fp)(...);
266 #else
267 typedef logical (*L_fp)();
268 #endif
269
270 static float spow_ui(float x, integer n) {
271         float pow=1.0; unsigned long int u;
272         if(n != 0) {
273                 if(n < 0) n = -n, x = 1/x;
274                 for(u = n; ; ) {
275                         if(u & 01) pow *= x;
276                         if(u >>= 1) x *= x;
277                         else break;
278                 }
279         }
280         return pow;
281 }
282 static double dpow_ui(double x, integer n) {
283         double pow=1.0; unsigned long int u;
284         if(n != 0) {
285                 if(n < 0) n = -n, x = 1/x;
286                 for(u = n; ; ) {
287                         if(u & 01) pow *= x;
288                         if(u >>= 1) x *= x;
289                         else break;
290                 }
291         }
292         return pow;
293 }
294 #ifdef _MSC_VER
295 static _Fcomplex cpow_ui(complex x, integer n) {
296         complex pow={1.0,0.0}; unsigned long int u;
297                 if(n != 0) {
298                 if(n < 0) n = -n, x.r = 1/x.r, x.i=1/x.i;
299                 for(u = n; ; ) {
300                         if(u & 01) pow.r *= x.r, pow.i *= x.i;
301                         if(u >>= 1) x.r *= x.r, x.i *= x.i;
302                         else break;
303                 }
304         }
305         _Fcomplex p={pow.r, pow.i};
306         return p;
307 }
308 #else
309 static _Complex float cpow_ui(_Complex float x, integer n) {
310         _Complex float pow=1.0; unsigned long int u;
311         if(n != 0) {
312                 if(n < 0) n = -n, x = 1/x;
313                 for(u = n; ; ) {
314                         if(u & 01) pow *= x;
315                         if(u >>= 1) x *= x;
316                         else break;
317                 }
318         }
319         return pow;
320 }
321 #endif
322 #ifdef _MSC_VER
323 static _Dcomplex zpow_ui(_Dcomplex x, integer n) {
324         _Dcomplex pow={1.0,0.0}; unsigned long int u;
325         if(n != 0) {
326                 if(n < 0) n = -n, x._Val[0] = 1/x._Val[0], x._Val[1] =1/x._Val[1];
327                 for(u = n; ; ) {
328                         if(u & 01) pow._Val[0] *= x._Val[0], pow._Val[1] *= x._Val[1];
329                         if(u >>= 1) x._Val[0] *= x._Val[0], x._Val[1] *= x._Val[1];
330                         else break;
331                 }
332         }
333         _Dcomplex p = {pow._Val[0], pow._Val[1]};
334         return p;
335 }
336 #else
337 static _Complex double zpow_ui(_Complex double x, integer n) {
338         _Complex double pow=1.0; unsigned long int u;
339         if(n != 0) {
340                 if(n < 0) n = -n, x = 1/x;
341                 for(u = n; ; ) {
342                         if(u & 01) pow *= x;
343                         if(u >>= 1) x *= x;
344                         else break;
345                 }
346         }
347         return pow;
348 }
349 #endif
350 static integer pow_ii(integer x, integer n) {
351         integer pow; unsigned long int u;
352         if (n <= 0) {
353                 if (n == 0 || x == 1) pow = 1;
354                 else if (x != -1) pow = x == 0 ? 1/x : 0;
355                 else n = -n;
356         }
357         if ((n > 0) || !(n == 0 || x == 1 || x != -1)) {
358                 u = n;
359                 for(pow = 1; ; ) {
360                         if(u & 01) pow *= x;
361                         if(u >>= 1) x *= x;
362                         else break;
363                 }
364         }
365         return pow;
366 }
367 static integer dmaxloc_(double *w, integer s, integer e, integer *n)
368 {
369         double m; integer i, mi;
370         for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
371                 if (w[i-1]>m) mi=i ,m=w[i-1];
372         return mi-s+1;
373 }
374 static integer smaxloc_(float *w, integer s, integer e, integer *n)
375 {
376         float m; integer i, mi;
377         for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
378                 if (w[i-1]>m) mi=i ,m=w[i-1];
379         return mi-s+1;
380 }
381 static inline void cdotc_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
382         integer n = *n_, incx = *incx_, incy = *incy_, i;
383 #ifdef _MSC_VER
384         _Fcomplex zdotc = {0.0, 0.0};
385         if (incx == 1 && incy == 1) {
386                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
387                         zdotc._Val[0] += conjf(Cf(&x[i]))._Val[0] * Cf(&y[i])._Val[0];
388                         zdotc._Val[1] += conjf(Cf(&x[i]))._Val[1] * Cf(&y[i])._Val[1];
389                 }
390         } else {
391                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
392                         zdotc._Val[0] += conjf(Cf(&x[i*incx]))._Val[0] * Cf(&y[i*incy])._Val[0];
393                         zdotc._Val[1] += conjf(Cf(&x[i*incx]))._Val[1] * Cf(&y[i*incy])._Val[1];
394                 }
395         }
396         pCf(z) = zdotc;
397 }
398 #else
399         _Complex float zdotc = 0.0;
400         if (incx == 1 && incy == 1) {
401                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
402                         zdotc += conjf(Cf(&x[i])) * Cf(&y[i]);
403                 }
404         } else {
405                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
406                         zdotc += conjf(Cf(&x[i*incx])) * Cf(&y[i*incy]);
407                 }
408         }
409         pCf(z) = zdotc;
410 }
411 #endif
412 static inline void zdotc_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
413         integer n = *n_, incx = *incx_, incy = *incy_, i;
414 #ifdef _MSC_VER
415         _Dcomplex zdotc = {0.0, 0.0};
416         if (incx == 1 && incy == 1) {
417                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
418                         zdotc._Val[0] += conj(Cd(&x[i]))._Val[0] * Cd(&y[i])._Val[0];
419                         zdotc._Val[1] += conj(Cd(&x[i]))._Val[1] * Cd(&y[i])._Val[1];
420                 }
421         } else {
422                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
423                         zdotc._Val[0] += conj(Cd(&x[i*incx]))._Val[0] * Cd(&y[i*incy])._Val[0];
424                         zdotc._Val[1] += conj(Cd(&x[i*incx]))._Val[1] * Cd(&y[i*incy])._Val[1];
425                 }
426         }
427         pCd(z) = zdotc;
428 }
429 #else
430         _Complex double zdotc = 0.0;
431         if (incx == 1 && incy == 1) {
432                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
433                         zdotc += conj(Cd(&x[i])) * Cd(&y[i]);
434                 }
435         } else {
436                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
437                         zdotc += conj(Cd(&x[i*incx])) * Cd(&y[i*incy]);
438                 }
439         }
440         pCd(z) = zdotc;
441 }
442 #endif  
443 static inline void cdotu_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
444         integer n = *n_, incx = *incx_, incy = *incy_, i;
445 #ifdef _MSC_VER
446         _Fcomplex zdotc = {0.0, 0.0};
447         if (incx == 1 && incy == 1) {
448                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
449                         zdotc._Val[0] += Cf(&x[i])._Val[0] * Cf(&y[i])._Val[0];
450                         zdotc._Val[1] += Cf(&x[i])._Val[1] * Cf(&y[i])._Val[1];
451                 }
452         } else {
453                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
454                         zdotc._Val[0] += Cf(&x[i*incx])._Val[0] * Cf(&y[i*incy])._Val[0];
455                         zdotc._Val[1] += Cf(&x[i*incx])._Val[1] * Cf(&y[i*incy])._Val[1];
456                 }
457         }
458         pCf(z) = zdotc;
459 }
460 #else
461         _Complex float zdotc = 0.0;
462         if (incx == 1 && incy == 1) {
463                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
464                         zdotc += Cf(&x[i]) * Cf(&y[i]);
465                 }
466         } else {
467                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
468                         zdotc += Cf(&x[i*incx]) * Cf(&y[i*incy]);
469                 }
470         }
471         pCf(z) = zdotc;
472 }
473 #endif
474 static inline void zdotu_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
475         integer n = *n_, incx = *incx_, incy = *incy_, i;
476 #ifdef _MSC_VER
477         _Dcomplex zdotc = {0.0, 0.0};
478         if (incx == 1 && incy == 1) {
479                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
480                         zdotc._Val[0] += Cd(&x[i])._Val[0] * Cd(&y[i])._Val[0];
481                         zdotc._Val[1] += Cd(&x[i])._Val[1] * Cd(&y[i])._Val[1];
482                 }
483         } else {
484                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
485                         zdotc._Val[0] += Cd(&x[i*incx])._Val[0] * Cd(&y[i*incy])._Val[0];
486                         zdotc._Val[1] += Cd(&x[i*incx])._Val[1] * Cd(&y[i*incy])._Val[1];
487                 }
488         }
489         pCd(z) = zdotc;
490 }
491 #else
492         _Complex double zdotc = 0.0;
493         if (incx == 1 && incy == 1) {
494                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
495                         zdotc += Cd(&x[i]) * Cd(&y[i]);
496                 }
497         } else {
498                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
499                         zdotc += Cd(&x[i*incx]) * Cd(&y[i*incy]);
500                 }
501         }
502         pCd(z) = zdotc;
503 }
504 #endif
505 /*  -- translated by f2c (version 20000121).
506    You must link the resulting object file with the libraries:
507         -lf2c -lm   (in that order)
508 */
509
510
511
512
513 /* Table of constant values */
514
515 static integer c_n1 = -1;
516 static logical c_false = FALSE_;
517
518 /* > \brief \b CUNCSD */
519
520 /*  =========== DOCUMENTATION =========== */
521
522 /* Online html documentation available at */
523 /*            http://www.netlib.org/lapack/explore-html/ */
524
525 /* > \htmlonly */
526 /* > Download CUNCSD + dependencies */
527 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cuncsd.
528 f"> */
529 /* > [TGZ]</a> */
530 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cuncsd.
531 f"> */
532 /* > [ZIP]</a> */
533 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cuncsd.
534 f"> */
535 /* > [TXT]</a> */
536 /* > \endhtmlonly */
537
538 /*  Definition: */
539 /*  =========== */
540
541 /*        SUBROUTINE CUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, */
542 /*                                    SIGNS, M, P, Q, X11, LDX11, X12, */
543 /*                                    LDX12, X21, LDX21, X22, LDX22, THETA, */
544 /*                                    U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, */
545 /*                                    LDV2T, WORK, LWORK, RWORK, LRWORK, */
546 /*                                    IWORK, INFO ) */
547
548 /*       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS */
549 /*       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12, */
550 /*      $                   LDX21, LDX22, LRWORK, LWORK, M, P, Q */
551 /*       INTEGER            IWORK( * ) */
552 /*       REAL               THETA( * ) */
553 /*       REAL               RWORK( * ) */
554 /*       COMPLEX            U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), */
555 /*      $                   V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ), */
556 /*      $                   X12( LDX12, * ), X21( LDX21, * ), X22( LDX22, */
557 /*      $                   * ) */
558
559
560 /* > \par Purpose: */
561 /*  ============= */
562 /* > */
563 /* > \verbatim */
564 /* > */
565 /* > CUNCSD computes the CS decomposition of an M-by-M partitioned */
566 /* > unitary matrix X: */
567 /* > */
568 /* >                                 [  I  0  0 |  0  0  0 ] */
569 /* >                                 [  0  C  0 |  0 -S  0 ] */
570 /* >     [ X11 | X12 ]   [ U1 |    ] [  0  0  0 |  0  0 -I ] [ V1 |    ]**H */
571 /* > X = [-----------] = [---------] [---------------------] [---------]   . */
572 /* >     [ X21 | X22 ]   [    | U2 ] [  0  0  0 |  I  0  0 ] [    | V2 ] */
573 /* >                                 [  0  S  0 |  0  C  0 ] */
574 /* >                                 [  0  0  I |  0  0  0 ] */
575 /* > */
576 /* > X11 is P-by-Q. The unitary matrices U1, U2, V1, and V2 are P-by-P, */
577 /* > (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are */
578 /* > R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in */
579 /* > which R = MIN(P,M-P,Q,M-Q). */
580 /* > \endverbatim */
581
582 /*  Arguments: */
583 /*  ========== */
584
585 /* > \param[in] JOBU1 */
586 /* > \verbatim */
587 /* >          JOBU1 is CHARACTER */
588 /* >          = 'Y':      U1 is computed; */
589 /* >          otherwise:  U1 is not computed. */
590 /* > \endverbatim */
591 /* > */
592 /* > \param[in] JOBU2 */
593 /* > \verbatim */
594 /* >          JOBU2 is CHARACTER */
595 /* >          = 'Y':      U2 is computed; */
596 /* >          otherwise:  U2 is not computed. */
597 /* > \endverbatim */
598 /* > */
599 /* > \param[in] JOBV1T */
600 /* > \verbatim */
601 /* >          JOBV1T is CHARACTER */
602 /* >          = 'Y':      V1T is computed; */
603 /* >          otherwise:  V1T is not computed. */
604 /* > \endverbatim */
605 /* > */
606 /* > \param[in] JOBV2T */
607 /* > \verbatim */
608 /* >          JOBV2T is CHARACTER */
609 /* >          = 'Y':      V2T is computed; */
610 /* >          otherwise:  V2T is not computed. */
611 /* > \endverbatim */
612 /* > */
613 /* > \param[in] TRANS */
614 /* > \verbatim */
615 /* >          TRANS is CHARACTER */
616 /* >          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major */
617 /* >                      order; */
618 /* >          otherwise:  X, U1, U2, V1T, and V2T are stored in column- */
619 /* >                      major order. */
620 /* > \endverbatim */
621 /* > */
622 /* > \param[in] SIGNS */
623 /* > \verbatim */
624 /* >          SIGNS is CHARACTER */
625 /* >          = 'O':      The lower-left block is made nonpositive (the */
626 /* >                      "other" convention); */
627 /* >          otherwise:  The upper-right block is made nonpositive (the */
628 /* >                      "default" convention). */
629 /* > \endverbatim */
630 /* > */
631 /* > \param[in] M */
632 /* > \verbatim */
633 /* >          M is INTEGER */
634 /* >          The number of rows and columns in X. */
635 /* > \endverbatim */
636 /* > */
637 /* > \param[in] P */
638 /* > \verbatim */
639 /* >          P is INTEGER */
640 /* >          The number of rows in X11 and X12. 0 <= P <= M. */
641 /* > \endverbatim */
642 /* > */
643 /* > \param[in] Q */
644 /* > \verbatim */
645 /* >          Q is INTEGER */
646 /* >          The number of columns in X11 and X21. 0 <= Q <= M. */
647 /* > \endverbatim */
648 /* > */
649 /* > \param[in,out] X11 */
650 /* > \verbatim */
651 /* >          X11 is COMPLEX array, dimension (LDX11,Q) */
652 /* >          On entry, part of the unitary matrix whose CSD is desired. */
653 /* > \endverbatim */
654 /* > */
655 /* > \param[in] LDX11 */
656 /* > \verbatim */
657 /* >          LDX11 is INTEGER */
658 /* >          The leading dimension of X11. LDX11 >= MAX(1,P). */
659 /* > \endverbatim */
660 /* > */
661 /* > \param[in,out] X12 */
662 /* > \verbatim */
663 /* >          X12 is COMPLEX array, dimension (LDX12,M-Q) */
664 /* >          On entry, part of the unitary matrix whose CSD is desired. */
665 /* > \endverbatim */
666 /* > */
667 /* > \param[in] LDX12 */
668 /* > \verbatim */
669 /* >          LDX12 is INTEGER */
670 /* >          The leading dimension of X12. LDX12 >= MAX(1,P). */
671 /* > \endverbatim */
672 /* > */
673 /* > \param[in,out] X21 */
674 /* > \verbatim */
675 /* >          X21 is COMPLEX array, dimension (LDX21,Q) */
676 /* >          On entry, part of the unitary matrix whose CSD is desired. */
677 /* > \endverbatim */
678 /* > */
679 /* > \param[in] LDX21 */
680 /* > \verbatim */
681 /* >          LDX21 is INTEGER */
682 /* >          The leading dimension of X11. LDX21 >= MAX(1,M-P). */
683 /* > \endverbatim */
684 /* > */
685 /* > \param[in,out] X22 */
686 /* > \verbatim */
687 /* >          X22 is COMPLEX array, dimension (LDX22,M-Q) */
688 /* >          On entry, part of the unitary matrix whose CSD is desired. */
689 /* > \endverbatim */
690 /* > */
691 /* > \param[in] LDX22 */
692 /* > \verbatim */
693 /* >          LDX22 is INTEGER */
694 /* >          The leading dimension of X11. LDX22 >= MAX(1,M-P). */
695 /* > \endverbatim */
696 /* > */
697 /* > \param[out] THETA */
698 /* > \verbatim */
699 /* >          THETA is REAL array, dimension (R), in which R = */
700 /* >          MIN(P,M-P,Q,M-Q). */
701 /* >          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and */
702 /* >          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ). */
703 /* > \endverbatim */
704 /* > */
705 /* > \param[out] U1 */
706 /* > \verbatim */
707 /* >          U1 is COMPLEX array, dimension (LDU1,P) */
708 /* >          If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1. */
709 /* > \endverbatim */
710 /* > */
711 /* > \param[in] LDU1 */
712 /* > \verbatim */
713 /* >          LDU1 is INTEGER */
714 /* >          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >= */
715 /* >          MAX(1,P). */
716 /* > \endverbatim */
717 /* > */
718 /* > \param[out] U2 */
719 /* > \verbatim */
720 /* >          U2 is COMPLEX array, dimension (LDU2,M-P) */
721 /* >          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary */
722 /* >          matrix U2. */
723 /* > \endverbatim */
724 /* > */
725 /* > \param[in] LDU2 */
726 /* > \verbatim */
727 /* >          LDU2 is INTEGER */
728 /* >          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >= */
729 /* >          MAX(1,M-P). */
730 /* > \endverbatim */
731 /* > */
732 /* > \param[out] V1T */
733 /* > \verbatim */
734 /* >          V1T is COMPLEX array, dimension (LDV1T,Q) */
735 /* >          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary */
736 /* >          matrix V1**H. */
737 /* > \endverbatim */
738 /* > */
739 /* > \param[in] LDV1T */
740 /* > \verbatim */
741 /* >          LDV1T is INTEGER */
742 /* >          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= */
743 /* >          MAX(1,Q). */
744 /* > \endverbatim */
745 /* > */
746 /* > \param[out] V2T */
747 /* > \verbatim */
748 /* >          V2T is COMPLEX array, dimension (LDV2T,M-Q) */
749 /* >          If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) unitary */
750 /* >          matrix V2**H. */
751 /* > \endverbatim */
752 /* > */
753 /* > \param[in] LDV2T */
754 /* > \verbatim */
755 /* >          LDV2T is INTEGER */
756 /* >          The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >= */
757 /* >          MAX(1,M-Q). */
758 /* > \endverbatim */
759 /* > */
760 /* > \param[out] WORK */
761 /* > \verbatim */
762 /* >          WORK is COMPLEX array, dimension (MAX(1,LWORK)) */
763 /* >          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
764 /* > \endverbatim */
765 /* > */
766 /* > \param[in] LWORK */
767 /* > \verbatim */
768 /* >          LWORK is INTEGER */
769 /* >          The dimension of the array WORK. */
770 /* > */
771 /* >          If LWORK = -1, then a workspace query is assumed; the routine */
772 /* >          only calculates the optimal size of the WORK array, returns */
773 /* >          this value as the first entry of the work array, and no error */
774 /* >          message related to LWORK is issued by XERBLA. */
775 /* > \endverbatim */
776 /* > */
777 /* > \param[out] RWORK */
778 /* > \verbatim */
779 /* >          RWORK is REAL array, dimension MAX(1,LRWORK) */
780 /* >          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. */
781 /* >          If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1), */
782 /* >          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R), */
783 /* >          define the matrix in intermediate bidiagonal-block form */
784 /* >          remaining after nonconvergence. INFO specifies the number */
785 /* >          of nonzero PHI's. */
786 /* > \endverbatim */
787 /* > */
788 /* > \param[in] LRWORK */
789 /* > \verbatim */
790 /* >          LRWORK is INTEGER */
791 /* >          The dimension of the array RWORK. */
792 /* > */
793 /* >          If LRWORK = -1, then a workspace query is assumed; the routine */
794 /* >          only calculates the optimal size of the RWORK array, returns */
795 /* >          this value as the first entry of the work array, and no error */
796 /* >          message related to LRWORK is issued by XERBLA. */
797 /* > \endverbatim */
798 /* > */
799 /* > \param[out] IWORK */
800 /* > \verbatim */
801 /* >          IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q)) */
802 /* > \endverbatim */
803 /* > */
804 /* > \param[out] INFO */
805 /* > \verbatim */
806 /* >          INFO is INTEGER */
807 /* >          = 0:  successful exit. */
808 /* >          < 0:  if INFO = -i, the i-th argument had an illegal value. */
809 /* >          > 0:  CBBCSD did not converge. See the description of RWORK */
810 /* >                above for details. */
811 /* > \endverbatim */
812
813 /* > \par References: */
814 /*  ================ */
815 /* > */
816 /* >  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. */
817 /* >      Algorithms, 50(1):33-65, 2009. */
818
819 /*  Authors: */
820 /*  ======== */
821
822 /* > \author Univ. of Tennessee */
823 /* > \author Univ. of California Berkeley */
824 /* > \author Univ. of Colorado Denver */
825 /* > \author NAG Ltd. */
826
827 /* > \date June 2016 */
828
829 /* > \ingroup complexOTHERcomputational */
830
831 /*  ===================================================================== */
832 /* Subroutine */ int cuncsd_(char *jobu1, char *jobu2, char *jobv1t, char *
833         jobv2t, char *trans, char *signs, integer *m, integer *p, integer *q, 
834         complex *x11, integer *ldx11, complex *x12, integer *ldx12, complex *
835         x21, integer *ldx21, complex *x22, integer *ldx22, real *theta, 
836         complex *u1, integer *ldu1, complex *u2, integer *ldu2, complex *v1t, 
837         integer *ldv1t, complex *v2t, integer *ldv2t, complex *work, integer *
838         lwork, real *rwork, integer *lrwork, integer *iwork, integer *info)
839 {
840     /* System generated locals */
841     integer u1_dim1, u1_offset, u2_dim1, u2_offset, v1t_dim1, v1t_offset, 
842             v2t_dim1, v2t_offset, x11_dim1, x11_offset, x12_dim1, x12_offset, 
843             x21_dim1, x21_offset, x22_dim1, x22_offset, i__1, i__2, i__3, 
844             i__4, i__5, i__6;
845
846     /* Local variables */
847     integer ib11d, ib11e, ib12d, ib12e, ib21d, ib21e, ib22d, ib22e, iphi;
848     logical colmajor;
849     integer lworkmin;
850     logical defaultsigns;
851     integer lworkopt, i__, j;
852     extern logical lsame_(char *, char *);
853     integer childinfo, p1, q1, lbbcsdworkmin, itaup1, itaup2, itauq1, itauq2, 
854             lorbdbworkmin, lrworkmin, lbbcsdworkopt;
855     logical wantu1, wantu2;
856     extern /* Subroutine */ int cbbcsd_(char *, char *, char *, char *, char *
857             , integer *, integer *, integer *, real *, real *, complex *, 
858             integer *, complex *, integer *, complex *, integer *, complex *, 
859             integer *, real *, real *, real *, real *, real *, real *, real *,
860              real *, real *, integer *, integer *);
861     integer lrworkopt, ibbcsd, lorbdbworkopt;
862     extern /* Subroutine */ int cunbdb_(char *, char *, integer *, integer *, 
863             integer *, complex *, integer *, complex *, integer *, complex *, 
864             integer *, complex *, integer *, real *, real *, complex *, 
865             complex *, complex *, complex *, complex *, integer *, integer *);
866     integer iorbdb, lorglqworkmin, lorgqrworkmin;
867     extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex 
868             *, integer *, complex *, integer *), xerbla_(char *, 
869             integer *, ftnlen), clapmr_(logical *, integer *, integer *, 
870             complex *, integer *, integer *), clapmt_(logical *, integer *, 
871             integer *, complex *, integer *, integer *);
872     integer lorglqworkopt;
873     extern /* Subroutine */ int cunglq_(integer *, integer *, integer *, 
874             complex *, integer *, complex *, complex *, integer *, integer *);
875     integer lorgqrworkopt, iorglq;
876     extern /* Subroutine */ int cungqr_(integer *, integer *, integer *, 
877             complex *, integer *, complex *, complex *, integer *, integer *);
878     integer iorgqr;
879     char signst[1], transt[1];
880     integer lbbcsdwork;
881     logical lquery;
882     integer lorbdbwork, lorglqwork, lorgqrwork;
883     logical wantv1t, wantv2t, lrquery;
884
885
886 /*  -- LAPACK computational routine (version 3.7.1) -- */
887 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
888 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
889 /*     June 2016 */
890
891
892 /*  =================================================================== */
893
894
895 /*     Test input arguments */
896
897     /* Parameter adjustments */
898     x11_dim1 = *ldx11;
899     x11_offset = 1 + x11_dim1 * 1;
900     x11 -= x11_offset;
901     x12_dim1 = *ldx12;
902     x12_offset = 1 + x12_dim1 * 1;
903     x12 -= x12_offset;
904     x21_dim1 = *ldx21;
905     x21_offset = 1 + x21_dim1 * 1;
906     x21 -= x21_offset;
907     x22_dim1 = *ldx22;
908     x22_offset = 1 + x22_dim1 * 1;
909     x22 -= x22_offset;
910     --theta;
911     u1_dim1 = *ldu1;
912     u1_offset = 1 + u1_dim1 * 1;
913     u1 -= u1_offset;
914     u2_dim1 = *ldu2;
915     u2_offset = 1 + u2_dim1 * 1;
916     u2 -= u2_offset;
917     v1t_dim1 = *ldv1t;
918     v1t_offset = 1 + v1t_dim1 * 1;
919     v1t -= v1t_offset;
920     v2t_dim1 = *ldv2t;
921     v2t_offset = 1 + v2t_dim1 * 1;
922     v2t -= v2t_offset;
923     --work;
924     --rwork;
925     --iwork;
926
927     /* Function Body */
928     *info = 0;
929     wantu1 = lsame_(jobu1, "Y");
930     wantu2 = lsame_(jobu2, "Y");
931     wantv1t = lsame_(jobv1t, "Y");
932     wantv2t = lsame_(jobv2t, "Y");
933     colmajor = ! lsame_(trans, "T");
934     defaultsigns = ! lsame_(signs, "O");
935     lquery = *lwork == -1;
936     lrquery = *lrwork == -1;
937     if (*m < 0) {
938         *info = -7;
939     } else if (*p < 0 || *p > *m) {
940         *info = -8;
941     } else if (*q < 0 || *q > *m) {
942         *info = -9;
943     } else if (colmajor && *ldx11 < f2cmax(1,*p)) {
944         *info = -11;
945     } else if (! colmajor && *ldx11 < f2cmax(1,*q)) {
946         *info = -11;
947     } else if (colmajor && *ldx12 < f2cmax(1,*p)) {
948         *info = -13;
949     } else /* if(complicated condition) */ {
950 /* Computing MAX */
951         i__1 = 1, i__2 = *m - *q;
952         if (! colmajor && *ldx12 < f2cmax(i__1,i__2)) {
953             *info = -13;
954         } else /* if(complicated condition) */ {
955 /* Computing MAX */
956             i__1 = 1, i__2 = *m - *p;
957             if (colmajor && *ldx21 < f2cmax(i__1,i__2)) {
958                 *info = -15;
959             } else if (! colmajor && *ldx21 < f2cmax(1,*q)) {
960                 *info = -15;
961             } else /* if(complicated condition) */ {
962 /* Computing MAX */
963                 i__1 = 1, i__2 = *m - *p;
964                 if (colmajor && *ldx22 < f2cmax(i__1,i__2)) {
965                     *info = -17;
966                 } else /* if(complicated condition) */ {
967 /* Computing MAX */
968                     i__1 = 1, i__2 = *m - *q;
969                     if (! colmajor && *ldx22 < f2cmax(i__1,i__2)) {
970                         *info = -17;
971                     } else if (wantu1 && *ldu1 < *p) {
972                         *info = -20;
973                     } else if (wantu2 && *ldu2 < *m - *p) {
974                         *info = -22;
975                     } else if (wantv1t && *ldv1t < *q) {
976                         *info = -24;
977                     } else if (wantv2t && *ldv2t < *m - *q) {
978                         *info = -26;
979                     }
980                 }
981             }
982         }
983     }
984
985 /*     Work with transpose if convenient */
986
987 /* Computing MIN */
988     i__1 = *p, i__2 = *m - *p;
989 /* Computing MIN */
990     i__3 = *q, i__4 = *m - *q;
991     if (*info == 0 && f2cmin(i__1,i__2) < f2cmin(i__3,i__4)) {
992         if (colmajor) {
993             *(unsigned char *)transt = 'T';
994         } else {
995             *(unsigned char *)transt = 'N';
996         }
997         if (defaultsigns) {
998             *(unsigned char *)signst = 'O';
999         } else {
1000             *(unsigned char *)signst = 'D';
1001         }
1002         cuncsd_(jobv1t, jobv2t, jobu1, jobu2, transt, signst, m, q, p, &x11[
1003                 x11_offset], ldx11, &x21[x21_offset], ldx21, &x12[x12_offset],
1004                  ldx12, &x22[x22_offset], ldx22, &theta[1], &v1t[v1t_offset], 
1005                 ldv1t, &v2t[v2t_offset], ldv2t, &u1[u1_offset], ldu1, &u2[
1006                 u2_offset], ldu2, &work[1], lwork, &rwork[1], lrwork, &iwork[
1007                 1], info);
1008         return 0;
1009     }
1010
1011 /*     Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if */
1012 /*     convenient */
1013
1014     if (*info == 0 && *m - *q < *q) {
1015         if (defaultsigns) {
1016             *(unsigned char *)signst = 'O';
1017         } else {
1018             *(unsigned char *)signst = 'D';
1019         }
1020         i__1 = *m - *p;
1021         i__2 = *m - *q;
1022         cuncsd_(jobu2, jobu1, jobv2t, jobv1t, trans, signst, m, &i__1, &i__2, 
1023                 &x22[x22_offset], ldx22, &x21[x21_offset], ldx21, &x12[
1024                 x12_offset], ldx12, &x11[x11_offset], ldx11, &theta[1], &u2[
1025                 u2_offset], ldu2, &u1[u1_offset], ldu1, &v2t[v2t_offset], 
1026                 ldv2t, &v1t[v1t_offset], ldv1t, &work[1], lwork, &rwork[1], 
1027                 lrwork, &iwork[1], info);
1028         return 0;
1029     }
1030
1031 /*     Compute workspace */
1032
1033     if (*info == 0) {
1034
1035 /*        Real workspace */
1036
1037         iphi = 2;
1038 /* Computing MAX */
1039         i__1 = 1, i__2 = *q - 1;
1040         ib11d = iphi + f2cmax(i__1,i__2);
1041         ib11e = ib11d + f2cmax(1,*q);
1042 /* Computing MAX */
1043         i__1 = 1, i__2 = *q - 1;
1044         ib12d = ib11e + f2cmax(i__1,i__2);
1045         ib12e = ib12d + f2cmax(1,*q);
1046 /* Computing MAX */
1047         i__1 = 1, i__2 = *q - 1;
1048         ib21d = ib12e + f2cmax(i__1,i__2);
1049         ib21e = ib21d + f2cmax(1,*q);
1050 /* Computing MAX */
1051         i__1 = 1, i__2 = *q - 1;
1052         ib22d = ib21e + f2cmax(i__1,i__2);
1053         ib22e = ib22d + f2cmax(1,*q);
1054 /* Computing MAX */
1055         i__1 = 1, i__2 = *q - 1;
1056         ibbcsd = ib22e + f2cmax(i__1,i__2);
1057         cbbcsd_(jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, &theta[1], &
1058                 theta[1], &u1[u1_offset], ldu1, &u2[u2_offset], ldu2, &v1t[
1059                 v1t_offset], ldv1t, &v2t[v2t_offset], ldv2t, &theta[1], &
1060                 theta[1], &theta[1], &theta[1], &theta[1], &theta[1], &theta[
1061                 1], &theta[1], &rwork[1], &c_n1, &childinfo);
1062         lbbcsdworkopt = (integer) rwork[1];
1063         lbbcsdworkmin = lbbcsdworkopt;
1064         lrworkopt = ibbcsd + lbbcsdworkopt - 1;
1065         lrworkmin = ibbcsd + lbbcsdworkmin - 1;
1066         rwork[1] = (real) lrworkopt;
1067
1068 /*        Complex workspace */
1069
1070         itaup1 = 2;
1071         itaup2 = itaup1 + f2cmax(1,*p);
1072 /* Computing MAX */
1073         i__1 = 1, i__2 = *m - *p;
1074         itauq1 = itaup2 + f2cmax(i__1,i__2);
1075         itauq2 = itauq1 + f2cmax(1,*q);
1076 /* Computing MAX */
1077         i__1 = 1, i__2 = *m - *q;
1078         iorgqr = itauq2 + f2cmax(i__1,i__2);
1079         i__1 = *m - *q;
1080         i__2 = *m - *q;
1081         i__3 = *m - *q;
1082 /* Computing MAX */
1083         i__5 = 1, i__6 = *m - *q;
1084         i__4 = f2cmax(i__5,i__6);
1085         cungqr_(&i__1, &i__2, &i__3, &u1[u1_offset], &i__4, &u1[u1_offset], &
1086                 work[1], &c_n1, &childinfo);
1087         lorgqrworkopt = (integer) work[1].r;
1088 /* Computing MAX */
1089         i__1 = 1, i__2 = *m - *q;
1090         lorgqrworkmin = f2cmax(i__1,i__2);
1091 /* Computing MAX */
1092         i__1 = 1, i__2 = *m - *q;
1093         iorglq = itauq2 + f2cmax(i__1,i__2);
1094         i__1 = *m - *q;
1095         i__2 = *m - *q;
1096         i__3 = *m - *q;
1097 /* Computing MAX */
1098         i__5 = 1, i__6 = *m - *q;
1099         i__4 = f2cmax(i__5,i__6);
1100         cunglq_(&i__1, &i__2, &i__3, &u1[u1_offset], &i__4, &u1[u1_offset], &
1101                 work[1], &c_n1, &childinfo);
1102         lorglqworkopt = (integer) work[1].r;
1103 /* Computing MAX */
1104         i__1 = 1, i__2 = *m - *q;
1105         lorglqworkmin = f2cmax(i__1,i__2);
1106 /* Computing MAX */
1107         i__1 = 1, i__2 = *m - *q;
1108         iorbdb = itauq2 + f2cmax(i__1,i__2);
1109         cunbdb_(trans, signs, m, p, q, &x11[x11_offset], ldx11, &x12[
1110                 x12_offset], ldx12, &x21[x21_offset], ldx21, &x22[x22_offset],
1111                  ldx22, &theta[1], &theta[1], &u1[u1_offset], &u2[u2_offset], 
1112                 &v1t[v1t_offset], &v2t[v2t_offset], &work[1], &c_n1, &
1113                 childinfo);
1114         lorbdbworkopt = (integer) work[1].r;
1115         lorbdbworkmin = lorbdbworkopt;
1116 /* Computing MAX */
1117         i__1 = iorgqr + lorgqrworkopt, i__2 = iorglq + lorglqworkopt, i__1 = 
1118                 f2cmax(i__1,i__2), i__2 = iorbdb + lorbdbworkopt;
1119         lworkopt = f2cmax(i__1,i__2) - 1;
1120 /* Computing MAX */
1121         i__1 = iorgqr + lorgqrworkmin, i__2 = iorglq + lorglqworkmin, i__1 = 
1122                 f2cmax(i__1,i__2), i__2 = iorbdb + lorbdbworkmin;
1123         lworkmin = f2cmax(i__1,i__2) - 1;
1124         i__1 = f2cmax(lworkopt,lworkmin);
1125         work[1].r = (real) i__1, work[1].i = 0.f;
1126
1127         if (*lwork < lworkmin && ! (lquery || lrquery)) {
1128             *info = -22;
1129         } else if (*lrwork < lrworkmin && ! (lquery || lrquery)) {
1130             *info = -24;
1131         } else {
1132             lorgqrwork = *lwork - iorgqr + 1;
1133             lorglqwork = *lwork - iorglq + 1;
1134             lorbdbwork = *lwork - iorbdb + 1;
1135             lbbcsdwork = *lrwork - ibbcsd + 1;
1136         }
1137     }
1138
1139 /*     Abort if any illegal arguments */
1140
1141     if (*info != 0) {
1142         i__1 = -(*info);
1143         xerbla_("CUNCSD", &i__1, (ftnlen)6);
1144         return 0;
1145     } else if (lquery || lrquery) {
1146         return 0;
1147     }
1148
1149 /*     Transform to bidiagonal block form */
1150
1151     cunbdb_(trans, signs, m, p, q, &x11[x11_offset], ldx11, &x12[x12_offset], 
1152             ldx12, &x21[x21_offset], ldx21, &x22[x22_offset], ldx22, &theta[1]
1153             , &rwork[iphi], &work[itaup1], &work[itaup2], &work[itauq1], &
1154             work[itauq2], &work[iorbdb], &lorbdbwork, &childinfo);
1155
1156 /*     Accumulate Householder reflectors */
1157
1158     if (colmajor) {
1159         if (wantu1 && *p > 0) {
1160             clacpy_("L", p, q, &x11[x11_offset], ldx11, &u1[u1_offset], ldu1);
1161             cungqr_(p, p, q, &u1[u1_offset], ldu1, &work[itaup1], &work[
1162                     iorgqr], &lorgqrwork, info);
1163         }
1164         if (wantu2 && *m - *p > 0) {
1165             i__1 = *m - *p;
1166             clacpy_("L", &i__1, q, &x21[x21_offset], ldx21, &u2[u2_offset], 
1167                     ldu2);
1168             i__1 = *m - *p;
1169             i__2 = *m - *p;
1170             cungqr_(&i__1, &i__2, q, &u2[u2_offset], ldu2, &work[itaup2], &
1171                     work[iorgqr], &lorgqrwork, info);
1172         }
1173         if (wantv1t && *q > 0) {
1174             i__1 = *q - 1;
1175             i__2 = *q - 1;
1176             clacpy_("U", &i__1, &i__2, &x11[(x11_dim1 << 1) + 1], ldx11, &v1t[
1177                     (v1t_dim1 << 1) + 2], ldv1t);
1178             i__1 = v1t_dim1 + 1;
1179             v1t[i__1].r = 1.f, v1t[i__1].i = 0.f;
1180             i__1 = *q;
1181             for (j = 2; j <= i__1; ++j) {
1182                 i__2 = j * v1t_dim1 + 1;
1183                 v1t[i__2].r = 0.f, v1t[i__2].i = 0.f;
1184                 i__2 = j + v1t_dim1;
1185                 v1t[i__2].r = 0.f, v1t[i__2].i = 0.f;
1186             }
1187             i__1 = *q - 1;
1188             i__2 = *q - 1;
1189             i__3 = *q - 1;
1190             cunglq_(&i__1, &i__2, &i__3, &v1t[(v1t_dim1 << 1) + 2], ldv1t, &
1191                     work[itauq1], &work[iorglq], &lorglqwork, info);
1192         }
1193         if (wantv2t && *m - *q > 0) {
1194             i__1 = *m - *q;
1195             clacpy_("U", p, &i__1, &x12[x12_offset], ldx12, &v2t[v2t_offset], 
1196                     ldv2t);
1197             if (*m - *p > *q) {
1198                 i__1 = *m - *p - *q;
1199                 i__2 = *m - *p - *q;
1200                 clacpy_("U", &i__1, &i__2, &x22[*q + 1 + (*p + 1) * x22_dim1],
1201                          ldx22, &v2t[*p + 1 + (*p + 1) * v2t_dim1], ldv2t);
1202             }
1203             if (*m > *q) {
1204                 i__1 = *m - *q;
1205                 i__2 = *m - *q;
1206                 i__3 = *m - *q;
1207                 cunglq_(&i__1, &i__2, &i__3, &v2t[v2t_offset], ldv2t, &work[
1208                         itauq2], &work[iorglq], &lorglqwork, info);
1209             }
1210         }
1211     } else {
1212         if (wantu1 && *p > 0) {
1213             clacpy_("U", q, p, &x11[x11_offset], ldx11, &u1[u1_offset], ldu1);
1214             cunglq_(p, p, q, &u1[u1_offset], ldu1, &work[itaup1], &work[
1215                     iorglq], &lorglqwork, info);
1216         }
1217         if (wantu2 && *m - *p > 0) {
1218             i__1 = *m - *p;
1219             clacpy_("U", q, &i__1, &x21[x21_offset], ldx21, &u2[u2_offset], 
1220                     ldu2);
1221             i__1 = *m - *p;
1222             i__2 = *m - *p;
1223             cunglq_(&i__1, &i__2, q, &u2[u2_offset], ldu2, &work[itaup2], &
1224                     work[iorglq], &lorglqwork, info);
1225         }
1226         if (wantv1t && *q > 0) {
1227             i__1 = *q - 1;
1228             i__2 = *q - 1;
1229             clacpy_("L", &i__1, &i__2, &x11[x11_dim1 + 2], ldx11, &v1t[(
1230                     v1t_dim1 << 1) + 2], ldv1t);
1231             i__1 = v1t_dim1 + 1;
1232             v1t[i__1].r = 1.f, v1t[i__1].i = 0.f;
1233             i__1 = *q;
1234             for (j = 2; j <= i__1; ++j) {
1235                 i__2 = j * v1t_dim1 + 1;
1236                 v1t[i__2].r = 0.f, v1t[i__2].i = 0.f;
1237                 i__2 = j + v1t_dim1;
1238                 v1t[i__2].r = 0.f, v1t[i__2].i = 0.f;
1239             }
1240             i__1 = *q - 1;
1241             i__2 = *q - 1;
1242             i__3 = *q - 1;
1243             cungqr_(&i__1, &i__2, &i__3, &v1t[(v1t_dim1 << 1) + 2], ldv1t, &
1244                     work[itauq1], &work[iorgqr], &lorgqrwork, info);
1245         }
1246         if (wantv2t && *m - *q > 0) {
1247 /* Computing MIN */
1248             i__1 = *p + 1;
1249             p1 = f2cmin(i__1,*m);
1250 /* Computing MIN */
1251             i__1 = *q + 1;
1252             q1 = f2cmin(i__1,*m);
1253             i__1 = *m - *q;
1254             clacpy_("L", &i__1, p, &x12[x12_offset], ldx12, &v2t[v2t_offset], 
1255                     ldv2t);
1256             if (*m > *p + *q) {
1257                 i__1 = *m - *p - *q;
1258                 i__2 = *m - *p - *q;
1259                 clacpy_("L", &i__1, &i__2, &x22[p1 + q1 * x22_dim1], ldx22, &
1260                         v2t[*p + 1 + (*p + 1) * v2t_dim1], ldv2t);
1261             }
1262             i__1 = *m - *q;
1263             i__2 = *m - *q;
1264             i__3 = *m - *q;
1265             cungqr_(&i__1, &i__2, &i__3, &v2t[v2t_offset], ldv2t, &work[
1266                     itauq2], &work[iorgqr], &lorgqrwork, info);
1267         }
1268     }
1269
1270 /*     Compute the CSD of the matrix in bidiagonal-block form */
1271
1272     cbbcsd_(jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, &theta[1], &rwork[
1273             iphi], &u1[u1_offset], ldu1, &u2[u2_offset], ldu2, &v1t[
1274             v1t_offset], ldv1t, &v2t[v2t_offset], ldv2t, &rwork[ib11d], &
1275             rwork[ib11e], &rwork[ib12d], &rwork[ib12e], &rwork[ib21d], &rwork[
1276             ib21e], &rwork[ib22d], &rwork[ib22e], &rwork[ibbcsd], &lbbcsdwork,
1277              info);
1278
1279 /*     Permute rows and columns to place identity submatrices in top- */
1280 /*     left corner of (1,1)-block and/or bottom-right corner of (1,2)- */
1281 /*     block and/or bottom-right corner of (2,1)-block and/or top-left */
1282 /*     corner of (2,2)-block */
1283
1284     if (*q > 0 && wantu2) {
1285         i__1 = *q;
1286         for (i__ = 1; i__ <= i__1; ++i__) {
1287             iwork[i__] = *m - *p - *q + i__;
1288         }
1289         i__1 = *m - *p;
1290         for (i__ = *q + 1; i__ <= i__1; ++i__) {
1291             iwork[i__] = i__ - *q;
1292         }
1293         if (colmajor) {
1294             i__1 = *m - *p;
1295             i__2 = *m - *p;
1296             clapmt_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
1297         } else {
1298             i__1 = *m - *p;
1299             i__2 = *m - *p;
1300             clapmr_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
1301         }
1302     }
1303     if (*m > 0 && wantv2t) {
1304         i__1 = *p;
1305         for (i__ = 1; i__ <= i__1; ++i__) {
1306             iwork[i__] = *m - *p - *q + i__;
1307         }
1308         i__1 = *m - *q;
1309         for (i__ = *p + 1; i__ <= i__1; ++i__) {
1310             iwork[i__] = i__ - *p;
1311         }
1312         if (! colmajor) {
1313             i__1 = *m - *q;
1314             i__2 = *m - *q;
1315             clapmt_(&c_false, &i__1, &i__2, &v2t[v2t_offset], ldv2t, &iwork[1]
1316                     );
1317         } else {
1318             i__1 = *m - *q;
1319             i__2 = *m - *q;
1320             clapmr_(&c_false, &i__1, &i__2, &v2t[v2t_offset], ldv2t, &iwork[1]
1321                     );
1322         }
1323     }
1324
1325     return 0;
1326
1327 /*     End CUNCSD */
1328
1329 } /* cuncsd_ */
1330