1*> \brief \b SORGTSQR
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SORGTSQR + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorgtsqr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorgtsqr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorgtsqr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SORGTSQR( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK,
22*      $                     INFO )
23*
24*       .. Scalar Arguments ..
25*       INTEGER           INFO, LDA, LDT, LWORK, M, N, MB, NB
26*       ..
27*       .. Array Arguments ..
28*       REAL              A( LDA, * ), T( LDT, * ), WORK( * )
29*       ..
30*
31*> \par Purpose:
32*  =============
33*>
34*> \verbatim
35*>
36*> SORGTSQR generates an M-by-N real matrix Q_out with orthonormal columns,
37*> which are the first N columns of a product of real orthogonal
38*> matrices of order M which are returned by SLATSQR
39*>
40*>      Q_out = first_N_columns_of( Q(1)_in * Q(2)_in * ... * Q(k)_in ).
41*>
42*> See the documentation for SLATSQR.
43*> \endverbatim
44*
45*  Arguments:
46*  ==========
47*
48*> \param[in] M
49*> \verbatim
50*>          M is INTEGER
51*>          The number of rows of the matrix A.  M >= 0.
52*> \endverbatim
53*>
54*> \param[in] N
55*> \verbatim
56*>          N is INTEGER
57*>          The number of columns of the matrix A. M >= N >= 0.
58*> \endverbatim
59*>
60*> \param[in] MB
61*> \verbatim
62*>          MB is INTEGER
63*>          The row block size used by SLATSQR to return
64*>          arrays A and T. MB > N.
65*>          (Note that if MB > M, then M is used instead of MB
66*>          as the row block size).
67*> \endverbatim
68*>
69*> \param[in] NB
70*> \verbatim
71*>          NB is INTEGER
72*>          The column block size used by SLATSQR to return
73*>          arrays A and T. NB >= 1.
74*>          (Note that if NB > N, then N is used instead of NB
75*>          as the column block size).
76*> \endverbatim
77*>
78*> \param[in,out] A
79*> \verbatim
80*>          A is REAL array, dimension (LDA,N)
81*>
82*>          On entry:
83*>
84*>             The elements on and above the diagonal are not accessed.
85*>             The elements below the diagonal represent the unit
86*>             lower-trapezoidal blocked matrix V computed by SLATSQR
87*>             that defines the input matrices Q_in(k) (ones on the
88*>             diagonal are not stored) (same format as the output A
89*>             below the diagonal in SLATSQR).
90*>
91*>          On exit:
92*>
93*>             The array A contains an M-by-N orthonormal matrix Q_out,
94*>             i.e the columns of A are orthogonal unit vectors.
95*> \endverbatim
96*>
97*> \param[in] LDA
98*> \verbatim
99*>          LDA is INTEGER
100*>          The leading dimension of the array A.  LDA >= max(1,M).
101*> \endverbatim
102*>
103*> \param[in] T
104*> \verbatim
105*>          T is REAL array,
106*>          dimension (LDT, N * NIRB)
107*>          where NIRB = Number_of_input_row_blocks
108*>                     = MAX( 1, CEIL((M-N)/(MB-N)) )
109*>          Let NICB = Number_of_input_col_blocks
110*>                   = CEIL(N/NB)
111*>
112*>          The upper-triangular block reflectors used to define the
113*>          input matrices Q_in(k), k=(1:NIRB*NICB). The block
114*>          reflectors are stored in compact form in NIRB block
115*>          reflector sequences. Each of NIRB block reflector sequences
116*>          is stored in a larger NB-by-N column block of T and consists
117*>          of NICB smaller NB-by-NB upper-triangular column blocks.
118*>          (same format as the output T in SLATSQR).
119*> \endverbatim
120*>
121*> \param[in] LDT
122*> \verbatim
123*>          LDT is INTEGER
124*>          The leading dimension of the array T.
125*>          LDT >= max(1,min(NB1,N)).
126*> \endverbatim
127*>
128*> \param[out] WORK
129*> \verbatim
130*>          (workspace) REAL array, dimension (MAX(2,LWORK))
131*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
132*> \endverbatim
133*>
134*> \param[in] LWORK
135*> \verbatim
136*>          The dimension of the array WORK.  LWORK >= (M+NB)*N.
137*>          If LWORK = -1, then a workspace query is assumed.
138*>          The routine only calculates the optimal size of the WORK
139*>          array, returns this value as the first entry of the WORK
140*>          array, and no error message related to LWORK is issued
141*>          by XERBLA.
142*> \endverbatim
143*>
144*> \param[out] INFO
145*> \verbatim
146*>          INFO is INTEGER
147*>          = 0:  successful exit
148*>          < 0:  if INFO = -i, the i-th argument had an illegal value
149*> \endverbatim
150*>
151*  Authors:
152*  ========
153*
154*> \author Univ. of Tennessee
155*> \author Univ. of California Berkeley
156*> \author Univ. of Colorado Denver
157*> \author NAG Ltd.
158*
159*> \ingroup singleOTHERcomputational
160*
161*> \par Contributors:
162*  ==================
163*>
164*> \verbatim
165*>
166*> November 2019, Igor Kozachenko,
167*>                Computer Science Division,
168*>                University of California, Berkeley
169*>
170*> \endverbatim
171*
172*  =====================================================================
173      SUBROUTINE SORGTSQR( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK,
174     $                     INFO )
175      IMPLICIT NONE
176*
177*  -- LAPACK computational routine --
178*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
179*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180*
181*     .. Scalar Arguments ..
182      INTEGER           INFO, LDA, LDT, LWORK, M, N, MB, NB
183*     ..
184*     .. Array Arguments ..
185      REAL              A( LDA, * ), T( LDT, * ), WORK( * )
186*     ..
187*
188*  =====================================================================
189*
190*     .. Parameters ..
191      REAL               ONE, ZERO
192      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
193*     ..
194*     .. Local Scalars ..
195      LOGICAL            LQUERY
196      INTEGER            IINFO, LDC, LWORKOPT, LC, LW, NBLOCAL, J
197*     ..
198*     .. External Subroutines ..
199      EXTERNAL           SCOPY, SLAMTSQR, SLASET, XERBLA
200*     ..
201*     .. Intrinsic Functions ..
202      INTRINSIC          REAL, MAX, MIN
203*     ..
204*     .. Executable Statements ..
205*
206*     Test the input parameters
207*
208      LQUERY  = LWORK.EQ.-1
209      INFO = 0
210      IF( M.LT.0 ) THEN
211         INFO = -1
212      ELSE IF( N.LT.0 .OR. M.LT.N ) THEN
213         INFO = -2
214      ELSE IF( MB.LE.N ) THEN
215         INFO = -3
216      ELSE IF( NB.LT.1 ) THEN
217         INFO = -4
218      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
219         INFO = -6
220      ELSE IF( LDT.LT.MAX( 1, MIN( NB, N ) ) ) THEN
221         INFO = -8
222      ELSE
223*
224*        Test the input LWORK for the dimension of the array WORK.
225*        This workspace is used to store array C(LDC, N) and WORK(LWORK)
226*        in the call to SLAMTSQR. See the documentation for SLAMTSQR.
227*
228         IF( LWORK.LT.2 .AND. (.NOT.LQUERY) ) THEN
229            INFO = -10
230         ELSE
231*
232*           Set block size for column blocks
233*
234            NBLOCAL = MIN( NB, N )
235*
236*           LWORK = -1, then set the size for the array C(LDC,N)
237*           in SLAMTSQR call and set the optimal size of the work array
238*           WORK(LWORK) in SLAMTSQR call.
239*
240            LDC = M
241            LC = LDC*N
242            LW = N * NBLOCAL
243*
244            LWORKOPT = LC+LW
245*
246            IF( ( LWORK.LT.MAX( 1, LWORKOPT ) ).AND.(.NOT.LQUERY) ) THEN
247               INFO = -10
248            END IF
249         END IF
250*
251      END IF
252*
253*     Handle error in the input parameters and return workspace query.
254*
255      IF( INFO.NE.0 ) THEN
256         CALL XERBLA( 'SORGTSQR', -INFO )
257         RETURN
258      ELSE IF ( LQUERY ) THEN
259         WORK( 1 ) = REAL( LWORKOPT )
260         RETURN
261      END IF
262*
263*     Quick return if possible
264*
265      IF( MIN( M, N ).EQ.0 ) THEN
266         WORK( 1 ) = REAL( LWORKOPT )
267         RETURN
268      END IF
269*
270*     (1) Form explicitly the tall-skinny M-by-N left submatrix Q1_in
271*     of M-by-M orthogonal matrix Q_in, which is implicitly stored in
272*     the subdiagonal part of input array A and in the input array T.
273*     Perform by the following operation using the routine SLAMTSQR.
274*
275*         Q1_in = Q_in * ( I ), where I is a N-by-N identity matrix,
276*                        ( 0 )        0 is a (M-N)-by-N zero matrix.
277*
278*     (1a) Form M-by-N matrix in the array WORK(1:LDC*N) with ones
279*     on the diagonal and zeros elsewhere.
280*
281      CALL SLASET( 'F', M, N, ZERO, ONE, WORK, LDC )
282*
283*     (1b)  On input, WORK(1:LDC*N) stores ( I );
284*                                          ( 0 )
285*
286*           On output, WORK(1:LDC*N) stores Q1_in.
287*
288      CALL SLAMTSQR( 'L', 'N', M, N, N, MB, NBLOCAL, A, LDA, T, LDT,
289     $               WORK, LDC, WORK( LC+1 ), LW, IINFO )
290*
291*     (2) Copy the result from the part of the work array (1:M,1:N)
292*     with the leading dimension LDC that starts at WORK(1) into
293*     the output array A(1:M,1:N) column-by-column.
294*
295      DO J = 1, N
296         CALL SCOPY( M, WORK( (J-1)*LDC + 1 ), 1, A( 1, J ), 1 )
297      END DO
298*
299      WORK( 1 ) = REAL( LWORKOPT )
300      RETURN
301*
302*     End of SORGTSQR
303*
304      END
305