LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
claqp2rk.f
Go to the documentation of this file.
1*> \brief \b CLAQP2RK computes truncated QR factorization with column pivoting of a complex matrix block using Level 2 BLAS and overwrites a complex m-by-nrhs matrix B with Q**H * B.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CLAQP2RK + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claqp2rk.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claqp2rk.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claqp2rk.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CLAQP2RK( M, N, NRHS, IOFFSET, KMAX, ABSTOL, RELTOL,
20* $ KP1, MAXC2NRM, A, LDA, K, MAXC2NRMK,
21* $ RELMAXC2NRMK, JPIV, TAU, VN1, VN2, WORK,
22* $ INFO )
23* IMPLICIT NONE
24*
25* .. Scalar Arguments ..
26* INTEGER INFO, IOFFSET, KP1, K, KMAX, LDA, M, N, NRHS
27* REAL ABSTOL, MAXC2NRM, MAXC2NRMK, RELMAXC2NRMK,
28* $ RELTOL
29* ..
30* .. Array Arguments ..
31* INTEGER JPIV( * )
32* REAL VN1( * ), VN2( * )
33* COMPLEX A( LDA, * ), TAU( * ), WORK( * )
34* $
35* ..
36*
37*
38*> \par Purpose:
39* =============
40*>
41*> \verbatim
42*>
43*> CLAQP2RK computes a truncated (rank K) or full rank Householder QR
44*> factorization with column pivoting of the complex matrix
45*> block A(IOFFSET+1:M,1:N) as
46*>
47*> A * P(K) = Q(K) * R(K).
48*>
49*> The routine uses Level 2 BLAS. The block A(1:IOFFSET,1:N)
50*> is accordingly pivoted, but not factorized.
51*>
52*> The routine also overwrites the right-hand-sides matrix block B
53*> stored in A(IOFFSET+1:M,N+1:N+NRHS) with Q(K)**H * B.
54*> \endverbatim
55*
56* Arguments:
57* ==========
58*
59*> \param[in] M
60*> \verbatim
61*> M is INTEGER
62*> The number of rows of the matrix A. M >= 0.
63*> \endverbatim
64*>
65*> \param[in] N
66*> \verbatim
67*> N is INTEGER
68*> The number of columns of the matrix A. N >= 0.
69*> \endverbatim
70*>
71*> \param[in] NRHS
72*> \verbatim
73*> NRHS is INTEGER
74*> The number of right hand sides, i.e., the number of
75*> columns of the matrix B. NRHS >= 0.
76*> \endverbatim
77*>
78*> \param[in] IOFFSET
79*> \verbatim
80*> IOFFSET is INTEGER
81*> The number of rows of the matrix A that must be pivoted
82*> but not factorized. IOFFSET >= 0.
83*>
84*> IOFFSET also represents the number of columns of the whole
85*> original matrix A_orig that have been factorized
86*> in the previous steps.
87*> \endverbatim
88*>
89*> \param[in] KMAX
90*> \verbatim
91*> KMAX is INTEGER
92*>
93*> The first factorization stopping criterion. KMAX >= 0.
94*>
95*> The maximum number of columns of the matrix A to factorize,
96*> i.e. the maximum factorization rank.
97*>
98*> a) If KMAX >= min(M-IOFFSET,N), then this stopping
99*> criterion is not used, factorize columns
100*> depending on ABSTOL and RELTOL.
101*>
102*> b) If KMAX = 0, then this stopping criterion is
103*> satisfied on input and the routine exits immediately.
104*> This means that the factorization is not performed,
105*> the matrices A and B and the arrays TAU, IPIV
106*> are not modified.
107*> \endverbatim
108*>
109*> \param[in] ABSTOL
110*> \verbatim
111*> ABSTOL is REAL, cannot be NaN.
112*>
113*> The second factorization stopping criterion.
114*>
115*> The absolute tolerance (stopping threshold) for
116*> maximum column 2-norm of the residual matrix.
117*> The algorithm converges (stops the factorization) when
118*> the maximum column 2-norm of the residual matrix
119*> is less than or equal to ABSTOL.
120*>
121*> a) If ABSTOL < 0.0, then this stopping criterion is not
122*> used, the routine factorizes columns depending
123*> on KMAX and RELTOL.
124*> This includes the case ABSTOL = -Inf.
125*>
126*> b) If 0.0 <= ABSTOL then the input value
127*> of ABSTOL is used.
128*> \endverbatim
129*>
130*> \param[in] RELTOL
131*> \verbatim
132*> RELTOL is REAL, cannot be NaN.
133*>
134*> The third factorization stopping criterion.
135*>
136*> The tolerance (stopping threshold) for the ratio of the
137*> maximum column 2-norm of the residual matrix to the maximum
138*> column 2-norm of the original matrix A_orig. The algorithm
139*> converges (stops the factorization), when this ratio is
140*> less than or equal to RELTOL.
141*>
142*> a) If RELTOL < 0.0, then this stopping criterion is not
143*> used, the routine factorizes columns depending
144*> on KMAX and ABSTOL.
145*> This includes the case RELTOL = -Inf.
146*>
147*> d) If 0.0 <= RELTOL then the input value of RELTOL
148*> is used.
149*> \endverbatim
150*>
151*> \param[in] KP1
152*> \verbatim
153*> KP1 is INTEGER
154*> The index of the column with the maximum 2-norm in
155*> the whole original matrix A_orig determined in the
156*> main routine CGEQP3RK. 1 <= KP1 <= N_orig_mat.
157*> \endverbatim
158*>
159*> \param[in] MAXC2NRM
160*> \verbatim
161*> MAXC2NRM is REAL
162*> The maximum column 2-norm of the whole original
163*> matrix A_orig computed in the main routine CGEQP3RK.
164*> MAXC2NRM >= 0.
165*> \endverbatim
166*>
167*> \param[in,out] A
168*> \verbatim
169*> A is COMPLEX array, dimension (LDA,N+NRHS)
170*> On entry:
171*> the M-by-N matrix A and M-by-NRHS matrix B, as in
172*>
173*> N NRHS
174*> array_A = M [ mat_A, mat_B ]
175*>
176*> On exit:
177*> 1. The elements in block A(IOFFSET+1:M,1:K) below
178*> the diagonal together with the array TAU represent
179*> the unitary matrix Q(K) as a product of elementary
180*> reflectors.
181*> 2. The upper triangular block of the matrix A stored
182*> in A(IOFFSET+1:M,1:K) is the triangular factor obtained.
183*> 3. The block of the matrix A stored in A(1:IOFFSET,1:N)
184*> has been accordingly pivoted, but not factorized.
185*> 4. The rest of the array A, block A(IOFFSET+1:M,K+1:N+NRHS).
186*> The left part A(IOFFSET+1:M,K+1:N) of this block
187*> contains the residual of the matrix A, and,
188*> if NRHS > 0, the right part of the block
189*> A(IOFFSET+1:M,N+1:N+NRHS) contains the block of
190*> the right-hand-side matrix B. Both these blocks have been
191*> updated by multiplication from the left by Q(K)**H.
192*> \endverbatim
193*>
194*> \param[in] LDA
195*> \verbatim
196*> LDA is INTEGER
197*> The leading dimension of the array A. LDA >= max(1,M).
198*> \endverbatim
199*>
200*> \param[out] K
201*> \verbatim
202*> K is INTEGER
203*> Factorization rank of the matrix A, i.e. the rank of
204*> the factor R, which is the same as the number of non-zero
205*> rows of the factor R. 0 <= K <= min(M-IOFFSET,KMAX,N).
206*>
207*> K also represents the number of non-zero Householder
208*> vectors.
209*> \endverbatim
210*>
211*> \param[out] MAXC2NRMK
212*> \verbatim
213*> MAXC2NRMK is REAL
214*> The maximum column 2-norm of the residual matrix,
215*> when the factorization stopped at rank K. MAXC2NRMK >= 0.
216*> \endverbatim
217*>
218*> \param[out] RELMAXC2NRMK
219*> \verbatim
220*> RELMAXC2NRMK is REAL
221*> The ratio MAXC2NRMK / MAXC2NRM of the maximum column
222*> 2-norm of the residual matrix (when the factorization
223*> stopped at rank K) to the maximum column 2-norm of the
224*> whole original matrix A. RELMAXC2NRMK >= 0.
225*> \endverbatim
226*>
227*> \param[out] JPIV
228*> \verbatim
229*> JPIV is INTEGER array, dimension (N)
230*> Column pivot indices, for 1 <= j <= N, column j
231*> of the matrix A was interchanged with column JPIV(j).
232*> \endverbatim
233*>
234*> \param[out] TAU
235*> \verbatim
236*> TAU is COMPLEX array, dimension (min(M-IOFFSET,N))
237*> The scalar factors of the elementary reflectors.
238*> \endverbatim
239*>
240*> \param[in,out] VN1
241*> \verbatim
242*> VN1 is REAL array, dimension (N)
243*> The vector with the partial column norms.
244*> \endverbatim
245*>
246*> \param[in,out] VN2
247*> \verbatim
248*> VN2 is REAL array, dimension (N)
249*> The vector with the exact column norms.
250*> \endverbatim
251*>
252*> \param[out] WORK
253*> \verbatim
254*> WORK is COMPLEX array, dimension (N-1)
255*> Used in CLARF1F subroutine to apply an elementary
256*> reflector from the left.
257*> \endverbatim
258*>
259*> \param[out] INFO
260*> \verbatim
261*> INFO is INTEGER
262*> 1) INFO = 0: successful exit.
263*> 2) If INFO = j_1, where 1 <= j_1 <= N, then NaN was
264*> detected and the routine stops the computation.
265*> The j_1-th column of the matrix A or the j_1-th
266*> element of array TAU contains the first occurrence
267*> of NaN in the factorization step K+1 ( when K columns
268*> have been factorized ).
269*>
270*> On exit:
271*> K is set to the number of
272*> factorized columns without
273*> exception.
274*> MAXC2NRMK is set to NaN.
275*> RELMAXC2NRMK is set to NaN.
276*> TAU(K+1:min(M,N)) is not set and contains undefined
277*> elements. If j_1=K+1, TAU(K+1)
278*> may contain NaN.
279*> 3) If INFO = j_2, where N+1 <= j_2 <= 2*N, then no NaN
280*> was detected, but +Inf (or -Inf) was detected and
281*> the routine continues the computation until completion.
282*> The (j_2-N)-th column of the matrix A contains the first
283*> occurrence of +Inf (or -Inf) in the factorization
284*> step K+1 ( when K columns have been factorized ).
285*> \endverbatim
286*
287* Authors:
288* ========
289*
290*> \author Univ. of Tennessee
291*> \author Univ. of California Berkeley
292*> \author Univ. of Colorado Denver
293*> \author NAG Ltd.
294*
295*> \ingroup laqp2rk
296*
297*> \par References:
298* ================
299*> [1] A Level 3 BLAS QR factorization algorithm with column pivoting developed in 1996.
300*> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain.
301*> X. Sun, Computer Science Dept., Duke University, USA.
302*> C. H. Bischof, Math. and Comp. Sci. Div., Argonne National Lab, USA.
303*> A BLAS-3 version of the QR factorization with column pivoting.
304*> LAPACK Working Note 114
305*> <a href="https://www.netlib.org/lapack/lawnspdf/lawn114.pdf">https://www.netlib.org/lapack/lawnspdf/lawn114.pdf</a>
306*> and in
307*> SIAM J. Sci. Comput., 19(5):1486-1494, Sept. 1998.
308*> <a href="https://doi.org/10.1137/S1064827595296732">https://doi.org/10.1137/S1064827595296732</a>
309*>
310*> [2] A partial column norm updating strategy developed in 2006.
311*> Z. Drmac and Z. Bujanovic, Dept. of Math., University of Zagreb, Croatia.
312*> On the failure of rank revealing QR factorization software – a case study.
313*> LAPACK Working Note 176.
314*> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">http://www.netlib.org/lapack/lawnspdf/lawn176.pdf</a>
315*> and in
316*> ACM Trans. Math. Softw. 35, 2, Article 12 (July 2008), 28 pages.
317*> <a href="https://doi.org/10.1145/1377612.1377616">https://doi.org/10.1145/1377612.1377616</a>
318*
319*> \par Contributors:
320* ==================
321*>
322*> \verbatim
323*>
324*> November 2023, Igor Kozachenko, James Demmel,
325*> EECS Department,
326*> University of California, Berkeley, USA.
327*>
328*> \endverbatim
329*
330* =====================================================================
331 SUBROUTINE claqp2rk( M, N, NRHS, IOFFSET, KMAX, ABSTOL, RELTOL,
332 $ KP1, MAXC2NRM, A, LDA, K, MAXC2NRMK,
333 $ RELMAXC2NRMK, JPIV, TAU, VN1, VN2, WORK,
334 $ INFO )
335 IMPLICIT NONE
336*
337* -- LAPACK auxiliary routine --
338* -- LAPACK is a software package provided by Univ. of Tennessee, --
339* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
340*
341* .. Scalar Arguments ..
342 INTEGER INFO, IOFFSET, KP1, K, KMAX, LDA, M, N, NRHS
343 REAL ABSTOL, MAXC2NRM, MAXC2NRMK, RELMAXC2NRMK,
344 $ RELTOL
345* ..
346* .. Array Arguments ..
347 INTEGER JPIV( * )
348 REAL VN1( * ), VN2( * )
349 COMPLEX A( LDA, * ), TAU( * ), WORK( * )
350* ..
351*
352* =====================================================================
353*
354* .. Parameters ..
355 REAL ZERO, ONE
356 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
357 COMPLEX CZERO
358 parameter( czero = ( 0.0e+0, 0.0e+0 ) )
359* ..
360* .. Local Scalars ..
361 INTEGER I, ITEMP, J, JMAXC2NRM, KK, KP, MINMNFACT,
362 $ MINMNUPDT
363 REAL HUGEVAL, TAUNAN, TEMP, TEMP2, TOL3Z
364* ..
365* .. External Subroutines ..
366 EXTERNAL clarf1f, clarfg, cswap
367* ..
368* .. Intrinsic Functions ..
369 INTRINSIC abs, real, conjg, aimag, max, min, sqrt
370* ..
371* .. External Functions ..
372 LOGICAL SISNAN
373 INTEGER ISAMAX
374 REAL SLAMCH, SCNRM2
375 EXTERNAL sisnan, slamch, isamax, scnrm2
376* ..
377* .. Executable Statements ..
378*
379* Initialize INFO
380*
381 info = 0
382*
383* MINMNFACT in the smallest dimension of the submatrix
384* A(IOFFSET+1:M,1:N) to be factorized.
385*
386* MINMNUPDT is the smallest dimension
387* of the subarray A(IOFFSET+1:M,1:N+NRHS) to be udated, which
388* contains the submatrices A(IOFFSET+1:M,1:N) and
389* B(IOFFSET+1:M,1:NRHS) as column blocks.
390*
391 minmnfact = min( m-ioffset, n )
392 minmnupdt = min( m-ioffset, n+nrhs )
393 kmax = min( kmax, minmnfact )
394 tol3z = sqrt( slamch( 'Epsilon' ) )
395 hugeval = slamch( 'Overflow' )
396*
397* Compute the factorization, KK is the lomn loop index.
398*
399 DO kk = 1, kmax
400*
401 i = ioffset + kk
402*
403 IF( i.EQ.1 ) THEN
404*
405* ============================================================
406*
407* We are at the first column of the original whole matrix A,
408* therefore we use the computed KP1 and MAXC2NRM from the
409* main routine.
410*
411 kp = kp1
412*
413* ============================================================
414*
415 ELSE
416*
417* ============================================================
418*
419* Determine the pivot column in KK-th step, i.e. the index
420* of the column with the maximum 2-norm in the
421* submatrix A(I:M,K:N).
422*
423 kp = ( kk-1 ) + isamax( n-kk+1, vn1( kk ), 1 )
424*
425* Determine the maximum column 2-norm and the relative maximum
426* column 2-norm of the submatrix A(I:M,KK:N) in step KK.
427* RELMAXC2NRMK will be computed later, after somecondition
428* checks on MAXC2NRMK.
429*
430 maxc2nrmk = vn1( kp )
431*
432* ============================================================
433*
434* Check if the submatrix A(I:M,KK:N) contains NaN, and set
435* INFO parameter to the column number, where the first NaN
436* is found and return from the routine.
437* We need to check the condition only if the
438* column index (same as row index) of the original whole
439* matrix is larger than 1, since the condition for whole
440* original matrix is checked in the main routine.
441*
442 IF( sisnan( maxc2nrmk ) ) THEN
443*
444* Set K, the number of factorized columns.
445* that are not zero.
446*
447 k = kk - 1
448 info = k + kp
449*
450* Set RELMAXC2NRMK to NaN.
451*
452 relmaxc2nrmk = maxc2nrmk
453*
454* Array TAU(K+1:MINMNFACT) is not set and contains
455* undefined elements.
456*
457 RETURN
458 END IF
459*
460* ============================================================
461*
462* Quick return, if the submatrix A(I:M,KK:N) is
463* a zero matrix.
464* We need to check the condition only if the
465* column index (same as row index) of the original whole
466* matrix is larger than 1, since the condition for whole
467* original matrix is checked in the main routine.
468*
469 IF( maxc2nrmk.EQ.zero ) THEN
470*
471* Set K, the number of factorized columns.
472* that are not zero.
473*
474 k = kk - 1
475 relmaxc2nrmk = zero
476*
477* Set TAUs corresponding to the columns that were not
478* factorized to ZERO, i.e. set TAU(KK:MINMNFACT) to CZERO.
479*
480 DO j = kk, minmnfact
481 tau( j ) = czero
482 END DO
483*
484* Return from the routine.
485*
486 RETURN
487*
488 END IF
489*
490* ============================================================
491*
492* Check if the submatrix A(I:M,KK:N) contains Inf,
493* set INFO parameter to the column number, where
494* the first Inf is found plus N, and continue
495* the computation.
496* We need to check the condition only if the
497* column index (same as row index) of the original whole
498* matrix is larger than 1, since the condition for whole
499* original matrix is checked in the main routine.
500*
501 IF( info.EQ.0 .AND. maxc2nrmk.GT.hugeval ) THEN
502 info = n + kk - 1 + kp
503 END IF
504*
505* ============================================================
506*
507* Test for the second and third stopping criteria.
508* NOTE: There is no need to test for ABSTOL >= ZERO, since
509* MAXC2NRMK is non-negative. Similarly, there is no need
510* to test for RELTOL >= ZERO, since RELMAXC2NRMK is
511* non-negative.
512* We need to check the condition only if the
513* column index (same as row index) of the original whole
514* matrix is larger than 1, since the condition for whole
515* original matrix is checked in the main routine.
516
517 relmaxc2nrmk = maxc2nrmk / maxc2nrm
518*
519 IF( maxc2nrmk.LE.abstol .OR. relmaxc2nrmk.LE.reltol ) THEN
520*
521* Set K, the number of factorized columns.
522*
523 k = kk - 1
524*
525* Set TAUs corresponding to the columns that were not
526* factorized to ZERO, i.e. set TAU(KK:MINMNFACT) to CZERO.
527*
528 DO j = kk, minmnfact
529 tau( j ) = czero
530 END DO
531*
532* Return from the routine.
533*
534 RETURN
535*
536 END IF
537*
538* ============================================================
539*
540* End ELSE of IF(I.EQ.1)
541*
542 END IF
543*
544* ===============================================================
545*
546* If the pivot column is not the first column of the
547* subblock A(1:M,KK:N):
548* 1) swap the KK-th column and the KP-th pivot column
549* in A(1:M,1:N);
550* 2) copy the KK-th element into the KP-th element of the partial
551* and exact 2-norm vectors VN1 and VN2. ( Swap is not needed
552* for VN1 and VN2 since we use the element with the index
553* larger than KK in the next loop step.)
554* 3) Save the pivot interchange with the indices relative to the
555* the original matrix A, not the block A(1:M,1:N).
556*
557 IF( kp.NE.kk ) THEN
558 CALL cswap( m, a( 1, kp ), 1, a( 1, kk ), 1 )
559 vn1( kp ) = vn1( kk )
560 vn2( kp ) = vn2( kk )
561 itemp = jpiv( kp )
562 jpiv( kp ) = jpiv( kk )
563 jpiv( kk ) = itemp
564 END IF
565*
566* Generate elementary reflector H(KK) using the column A(I:M,KK),
567* if the column has more than one element, otherwise
568* the elementary reflector would be an identity matrix,
569* and TAU(KK) = CZERO.
570*
571 IF( i.LT.m ) THEN
572 CALL clarfg( m-i+1, a( i, kk ), a( i+1, kk ), 1,
573 $ tau( kk ) )
574 ELSE
575 tau( kk ) = czero
576 END IF
577*
578* Check if TAU(KK) contains NaN, set INFO parameter
579* to the column number where NaN is found and return from
580* the routine.
581* NOTE: There is no need to check TAU(KK) for Inf,
582* since CLARFG cannot produce TAU(KK) or Householder vector
583* below the diagonal containing Inf. Only BETA on the diagonal,
584* returned by CLARFG can contain Inf, which requires
585* TAU(KK) to contain NaN. Therefore, this case of generating Inf
586* by CLARFG is covered by checking TAU(KK) for NaN.
587*
588 IF( sisnan( real( tau(kk) ) ) ) THEN
589 taunan = real( tau(kk) )
590 ELSE IF( sisnan( aimag( tau(kk) ) ) ) THEN
591 taunan = aimag( tau(kk) )
592 ELSE
593 taunan = zero
594 END IF
595*
596 IF( sisnan( taunan ) ) THEN
597 k = kk - 1
598 info = kk
599*
600* Set MAXC2NRMK and RELMAXC2NRMK to NaN.
601*
602 maxc2nrmk = taunan
603 relmaxc2nrmk = taunan
604*
605* Array TAU(KK:MINMNFACT) is not set and contains
606* undefined elements, except the first element TAU(KK) = NaN.
607*
608 RETURN
609 END IF
610*
611* Apply H(KK)**H to A(I:M,KK+1:N+NRHS) from the left.
612* ( If M >= N, then at KK = N there is no residual matrix,
613* i.e. no columns of A to update, only columns of B.
614* If M < N, then at KK = M-IOFFSET, I = M and we have a
615* one-row residual matrix in A and the elementary
616* reflector is a unit matrix, TAU(KK) = CZERO, i.e. no update
617* is needed for the residual matrix in A and the
618* right-hand-side-matrix in B.
619* Therefore, we update only if
620* KK < MINMNUPDT = min(M-IOFFSET, N+NRHS)
621* condition is satisfied, not only KK < N+NRHS )
622*
623 IF( kk.LT.minmnupdt ) THEN
624 CALL clarf1f( 'Left', m-i+1, n+nrhs-kk, a( i, kk ), 1,
625 $ conjg( tau( kk ) ), a( i, kk+1 ), lda,
626 $ work( 1 ) )
627 END IF
628*
629 IF( kk.LT.minmnfact ) THEN
630*
631* Update the partial column 2-norms for the residual matrix,
632* only if the residual matrix A(I+1:M,KK+1:N) exists, i.e.
633* when KK < min(M-IOFFSET, N).
634*
635 DO j = kk + 1, n
636 IF( vn1( j ).NE.zero ) THEN
637*
638* NOTE: The following lines follow from the analysis in
639* Lapack Working Note 176.
640*
641 temp = one - ( abs( a( i, j ) ) / vn1( j ) )**2
642 temp = max( temp, zero )
643 temp2 = temp*( vn1( j ) / vn2( j ) )**2
644 IF( temp2 .LE. tol3z ) THEN
645*
646* Compute the column 2-norm for the partial
647* column A(I+1:M,J) by explicitly computing it,
648* and store it in both partial 2-norm vector VN1
649* and exact column 2-norm vector VN2.
650*
651 vn1( j ) = scnrm2( m-i, a( i+1, j ), 1 )
652 vn2( j ) = vn1( j )
653*
654 ELSE
655*
656* Update the column 2-norm for the partial
657* column A(I+1:M,J) by removing one
658* element A(I,J) and store it in partial
659* 2-norm vector VN1.
660*
661 vn1( j ) = vn1( j )*sqrt( temp )
662*
663 END IF
664 END IF
665 END DO
666*
667 END IF
668*
669* End factorization loop
670*
671 END DO
672*
673* If we reached this point, all colunms have been factorized,
674* i.e. no condition was triggered to exit the routine.
675* Set the number of factorized columns.
676*
677 k = kmax
678*
679* We reached the end of the loop, i.e. all KMAX columns were
680* factorized, we need to set MAXC2NRMK and RELMAXC2NRMK before
681* we return.
682*
683 IF( k.LT.minmnfact ) THEN
684*
685 jmaxc2nrm = k + isamax( n-k, vn1( k+1 ), 1 )
686 maxc2nrmk = vn1( jmaxc2nrm )
687*
688 IF( k.EQ.0 ) THEN
689 relmaxc2nrmk = one
690 ELSE
691 relmaxc2nrmk = maxc2nrmk / maxc2nrm
692 END IF
693*
694 ELSE
695 maxc2nrmk = zero
696 relmaxc2nrmk = zero
697 END IF
698*
699* We reached the end of the loop, i.e. all KMAX columns were
700* factorized, set TAUs corresponding to the columns that were
701* not factorized to ZERO, i.e. TAU(K+1:MINMNFACT) set to CZERO.
702*
703 DO j = k + 1, minmnfact
704 tau( j ) = czero
705 END DO
706*
707 RETURN
708*
709* End of CLAQP2RK
710*
711 END
subroutine claqp2rk(m, n, nrhs, ioffset, kmax, abstol, reltol, kp1, maxc2nrm, a, lda, k, maxc2nrmk, relmaxc2nrmk, jpiv, tau, vn1, vn2, work, info)
CLAQP2RK computes truncated QR factorization with column pivoting of a complex matrix block using Lev...
Definition claqp2rk.f:335
subroutine clarf1f(side, m, n, v, incv, tau, c, ldc, work)
CLARF1F applies an elementary reflector to a general rectangular
Definition clarf1f.f:126
subroutine clarfg(n, alpha, x, incx, tau)
CLARFG generates an elementary reflector (Householder matrix).
Definition clarfg.f:104
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81