LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
claqp3rk.f
Go to the documentation of this file.
1*> \brief \b CLAQP3RK computes a step of truncated QR factorization with column pivoting of a complex m-by-n matrix A using Level 3 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*> \htmlonly
9*> Download CLAQP3RK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claqp3rk.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claqp3rk.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claqp3rk.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLAQP3RK( M, N, NRHS, IOFFSET, NB, ABSTOL,
22* $ RELTOL, KP1, MAXC2NRM, A, LDA, DONE, KB,
23* $ MAXC2NRMK, RELMAXC2NRMK, JPIV, TAU,
24* $ VN1, VN2, AUXV, F, LDF, IWORK, INFO )
25* IMPLICIT NONE
26* LOGICAL DONE
27* INTEGER INFO, IOFFSET, KB, KP1, LDA, LDF, M, N,
28* $ NB, NRHS
29* REAL ABSTOL, MAXC2NRM, MAXC2NRMK, RELMAXC2NRMK,
30* $ RELTOL
31* ..
32* .. Array Arguments ..
33* INTEGER IWORK( * ), JPIV( * )
34* REAL VN1( * ), VN2( * )
35* COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * )
36* ..
37*
38*
39*> \par Purpose:
40* =============
41*>
42*> \verbatim
43*>
44*> CLAQP3RK computes a step of truncated QR factorization with column
45*> pivoting of a complex M-by-N matrix A block A(IOFFSET+1:M,1:N)
46*> by using Level 3 BLAS as
47*>
48*> A * P(KB) = Q(KB) * R(KB).
49*>
50*> The routine tries to factorize NB columns from A starting from
51*> the row IOFFSET+1 and updates the residual matrix with BLAS 3
52*> xGEMM. The number of actually factorized columns is returned
53*> is smaller than NB.
54*>
55*> Block A(1:IOFFSET,1:N) is accordingly pivoted, but not factorized.
56*>
57*> The routine also overwrites the right-hand-sides B matrix stored
58*> in A(IOFFSET+1:M,1:N+1:N+NRHS) with Q(KB)**H * B.
59*>
60*> Cases when the number of factorized columns KB < NB:
61*>
62*> (1) In some cases, due to catastrophic cancellations, it cannot
63*> factorize all NB columns and need to update the residual matrix.
64*> Hence, the actual number of factorized columns in the block returned
65*> in KB is smaller than NB. The logical DONE is returned as FALSE.
66*> The factorization of the whole original matrix A_orig must proceed
67*> with the next block.
68*>
69*> (2) Whenever the stopping criterion ABSTOL or RELTOL is satisfied,
70*> the factorization of the whole original matrix A_orig is stopped,
71*> the logical DONE is returned as TRUE. The number of factorized
72*> columns which is smaller than NB is returned in KB.
73*>
74*> (3) In case both stopping criteria ABSTOL or RELTOL are not used,
75*> and when the residual matrix is a zero matrix in some factorization
76*> step KB, the factorization of the whole original matrix A_orig is
77*> stopped, the logical DONE is returned as TRUE. The number of
78*> factorized columns which is smaller than NB is returned in KB.
79*>
80*> (4) Whenever NaN is detected in the matrix A or in the array TAU,
81*> the factorization of the whole original matrix A_orig is stopped,
82*> the logical DONE is returned as TRUE. The number of factorized
83*> columns which is smaller than NB is returned in KB. The INFO
84*> parameter is set to the column index of the first NaN occurrence.
85*>
86*> \endverbatim
87*
88* Arguments:
89* ==========
90*
91*> \param[in] M
92*> \verbatim
93*> M is INTEGER
94*> The number of rows of the matrix A. M >= 0.
95*> \endverbatim
96*>
97*> \param[in] N
98*> \verbatim
99*> N is INTEGER
100*> The number of columns of the matrix A. N >= 0
101*> \endverbatim
102*>
103*> \param[in] NRHS
104*> \verbatim
105*> NRHS is INTEGER
106*> The number of right hand sides, i.e., the number of
107*> columns of the matrix B. NRHS >= 0.
108*> \endverbatim
109*>
110*> \param[in] IOFFSET
111*> \verbatim
112*> IOFFSET is INTEGER
113*> The number of rows of the matrix A that must be pivoted
114*> but not factorized. IOFFSET >= 0.
115*>
116*> IOFFSET also represents the number of columns of the whole
117*> original matrix A_orig that have been factorized
118*> in the previous steps.
119*> \endverbatim
120*>
121*> \param[in] NB
122*> \verbatim
123*> NB is INTEGER
124*> Factorization block size, i.e the number of columns
125*> to factorize in the matrix A. 0 <= NB
126*>
127*> If NB = 0, then the routine exits immediately.
128*> This means that the factorization is not performed,
129*> the matrices A and B and the arrays TAU, IPIV
130*> are not modified.
131*> \endverbatim
132*>
133*> \param[in] ABSTOL
134*> \verbatim
135*> ABSTOL is REAL, cannot be NaN.
136*>
137*> The absolute tolerance (stopping threshold) for
138*> maximum column 2-norm of the residual matrix.
139*> The algorithm converges (stops the factorization) when
140*> the maximum column 2-norm of the residual matrix
141*> is less than or equal to ABSTOL.
142*>
143*> a) If ABSTOL < 0.0, then this stopping criterion is not
144*> used, the routine factorizes columns depending
145*> on NB and RELTOL.
146*> This includes the case ABSTOL = -Inf.
147*>
148*> b) If 0.0 <= ABSTOL then the input value
149*> of ABSTOL is used.
150*> \endverbatim
151*>
152*> \param[in] RELTOL
153*> \verbatim
154*> RELTOL is REAL, cannot be NaN.
155*>
156*> The tolerance (stopping threshold) for the ratio of the
157*> maximum column 2-norm of the residual matrix to the maximum
158*> column 2-norm of the original matrix A_orig. The algorithm
159*> converges (stops the factorization), when this ratio is
160*> less than or equal to RELTOL.
161*>
162*> a) If RELTOL < 0.0, then this stopping criterion is not
163*> used, the routine factorizes columns depending
164*> on NB and ABSTOL.
165*> This includes the case RELTOL = -Inf.
166*>
167*> d) If 0.0 <= RELTOL then the input value of RELTOL
168*> is used.
169*> \endverbatim
170*>
171*> \param[in] KP1
172*> \verbatim
173*> KP1 is INTEGER
174*> The index of the column with the maximum 2-norm in
175*> the whole original matrix A_orig determined in the
176*> main routine CGEQP3RK. 1 <= KP1 <= N_orig.
177*> \endverbatim
178*>
179*> \param[in] MAXC2NRM
180*> \verbatim
181*> MAXC2NRM is REAL
182*> The maximum column 2-norm of the whole original
183*> matrix A_orig computed in the main routine CGEQP3RK.
184*> MAXC2NRM >= 0.
185*> \endverbatim
186*>
187*> \param[in,out] A
188*> \verbatim
189*> A is COMPLEX array, dimension (LDA,N+NRHS)
190*> On entry:
191*> the M-by-N matrix A and M-by-NRHS matrix B, as in
192*>
193*> N NRHS
194*> array_A = M [ mat_A, mat_B ]
195*>
196*> On exit:
197*> 1. The elements in block A(IOFFSET+1:M,1:KB) below
198*> the diagonal together with the array TAU represent
199*> the unitary matrix Q(KB) as a product of elementary
200*> reflectors.
201*> 2. The upper triangular block of the matrix A stored
202*> in A(IOFFSET+1:M,1:KB) is the triangular factor obtained.
203*> 3. The block of the matrix A stored in A(1:IOFFSET,1:N)
204*> has been accordingly pivoted, but not factorized.
205*> 4. The rest of the array A, block A(IOFFSET+1:M,KB+1:N+NRHS).
206*> The left part A(IOFFSET+1:M,KB+1:N) of this block
207*> contains the residual of the matrix A, and,
208*> if NRHS > 0, the right part of the block
209*> A(IOFFSET+1:M,N+1:N+NRHS) contains the block of
210*> the right-hand-side matrix B. Both these blocks have been
211*> updated by multiplication from the left by Q(KB)**H.
212*> \endverbatim
213*>
214*> \param[in] LDA
215*> \verbatim
216*> LDA is INTEGER
217*> The leading dimension of the array A. LDA >= max(1,M).
218*> \endverbatim
219*>
220*> \param[out]
221*> \verbatim
222*> DONE is LOGICAL
223*> TRUE: a) if the factorization completed before processing
224*> all min(M-IOFFSET,NB,N) columns due to ABSTOL
225*> or RELTOL criterion,
226*> b) if the factorization completed before processing
227*> all min(M-IOFFSET,NB,N) columns due to the
228*> residual matrix being a ZERO matrix.
229*> c) when NaN was detected in the matrix A
230*> or in the array TAU.
231*> FALSE: otherwise.
232*> \endverbatim
233*>
234*> \param[out] KB
235*> \verbatim
236*> KB is INTEGER
237*> Factorization rank of the matrix A, i.e. the rank of
238*> the factor R, which is the same as the number of non-zero
239*> rows of the factor R. 0 <= KB <= min(M-IOFFSET,NB,N).
240*>
241*> KB also represents the number of non-zero Householder
242*> vectors.
243*> \endverbatim
244*>
245*> \param[out] MAXC2NRMK
246*> \verbatim
247*> MAXC2NRMK is REAL
248*> The maximum column 2-norm of the residual matrix,
249*> when the factorization stopped at rank KB. MAXC2NRMK >= 0.
250*> \endverbatim
251*>
252*> \param[out] RELMAXC2NRMK
253*> \verbatim
254*> RELMAXC2NRMK is REAL
255*> The ratio MAXC2NRMK / MAXC2NRM of the maximum column
256*> 2-norm of the residual matrix (when the factorization
257*> stopped at rank KB) to the maximum column 2-norm of the
258*> original matrix A_orig. RELMAXC2NRMK >= 0.
259*> \endverbatim
260*>
261*> \param[out] JPIV
262*> \verbatim
263*> JPIV is INTEGER array, dimension (N)
264*> Column pivot indices, for 1 <= j <= N, column j
265*> of the matrix A was interchanged with column JPIV(j).
266*> \endverbatim
267*>
268*> \param[out] TAU
269*> \verbatim
270*> TAU is COMPLEX array, dimension (min(M-IOFFSET,N))
271*> The scalar factors of the elementary reflectors.
272*> \endverbatim
273*>
274*> \param[in,out] VN1
275*> \verbatim
276*> VN1 is REAL array, dimension (N)
277*> The vector with the partial column norms.
278*> \endverbatim
279*>
280*> \param[in,out] VN2
281*> \verbatim
282*> VN2 is REAL array, dimension (N)
283*> The vector with the exact column norms.
284*> \endverbatim
285*>
286*> \param[out] AUXV
287*> \verbatim
288*> AUXV is COMPLEX array, dimension (NB)
289*> Auxiliary vector.
290*> \endverbatim
291*>
292*> \param[out] F
293*> \verbatim
294*> F is COMPLEX array, dimension (LDF,NB)
295*> Matrix F**H = L*(Y**H)*A.
296*> \endverbatim
297*>
298*> \param[in] LDF
299*> \verbatim
300*> LDF is INTEGER
301*> The leading dimension of the array F. LDF >= max(1,N+NRHS).
302*> \endverbatim
303*>
304*> \param[out] IWORK
305*> \verbatim
306*> IWORK is INTEGER array, dimension (N-1).
307*> Is a work array. ( IWORK is used to store indices
308*> of "bad" columns for norm downdating in the residual
309*> matrix ).
310*> \endverbatim
311*>
312*> \param[out] INFO
313*> \verbatim
314*> INFO is INTEGER
315*> 1) INFO = 0: successful exit.
316*> 2) If INFO = j_1, where 1 <= j_1 <= N, then NaN was
317*> detected and the routine stops the computation.
318*> The j_1-th column of the matrix A or the j_1-th
319*> element of array TAU contains the first occurrence
320*> of NaN in the factorization step KB+1 ( when KB columns
321*> have been factorized ).
322*>
323*> On exit:
324*> KB is set to the number of
325*> factorized columns without
326*> exception.
327*> MAXC2NRMK is set to NaN.
328*> RELMAXC2NRMK is set to NaN.
329*> TAU(KB+1:min(M,N)) is not set and contains undefined
330*> elements. If j_1=KB+1, TAU(KB+1)
331*> may contain NaN.
332*> 3) If INFO = j_2, where N+1 <= j_2 <= 2*N, then no NaN
333*> was detected, but +Inf (or -Inf) was detected and
334*> the routine continues the computation until completion.
335*> The (j_2-N)-th column of the matrix A contains the first
336*> occurrence of +Inf (or -Inf) in the actorization
337*> step KB+1 ( when KB columns have been factorized ).
338*> \endverbatim
339*
340* Authors:
341* ========
342*
343*> \author Univ. of Tennessee
344*> \author Univ. of California Berkeley
345*> \author Univ. of Colorado Denver
346*> \author NAG Ltd.
347*
348*> \ingroup laqp3rk
349*
350*> \par References:
351* ================
352*> [1] A Level 3 BLAS QR factorization algorithm with column pivoting developed in 1996.
353*> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain.
354*> X. Sun, Computer Science Dept., Duke University, USA.
355*> C. H. Bischof, Math. and Comp. Sci. Div., Argonne National Lab, USA.
356*> A BLAS-3 version of the QR factorization with column pivoting.
357*> LAPACK Working Note 114
358*> \htmlonly
359*> <a href="https://www.netlib.org/lapack/lawnspdf/lawn114.pdf">https://www.netlib.org/lapack/lawnspdf/lawn114.pdf</a>
360*> \endhtmlonly
361*> and in
362*> SIAM J. Sci. Comput., 19(5):1486-1494, Sept. 1998.
363*> \htmlonly
364*> <a href="https://doi.org/10.1137/S1064827595296732">https://doi.org/10.1137/S1064827595296732</a>
365*> \endhtmlonly
366*>
367*> [2] A partial column norm updating strategy developed in 2006.
368*> Z. Drmac and Z. Bujanovic, Dept. of Math., University of Zagreb, Croatia.
369*> On the failure of rank revealing QR factorization software – a case study.
370*> LAPACK Working Note 176.
371*> \htmlonly
372*> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">http://www.netlib.org/lapack/lawnspdf/lawn176.pdf</a>
373*> \endhtmlonly
374*> and in
375*> ACM Trans. Math. Softw. 35, 2, Article 12 (July 2008), 28 pages.
376*> \htmlonly
377*> <a href="https://doi.org/10.1145/1377612.1377616">https://doi.org/10.1145/1377612.1377616</a>
378*> \endhtmlonly
379*
380*> \par Contributors:
381* ==================
382*>
383*> \verbatim
384*>
385*> November 2023, Igor Kozachenko, James Demmel,
386*> EECS Department,
387*> University of California, Berkeley, USA.
388*>
389*> \endverbatim
390*
391* =====================================================================
392 SUBROUTINE claqp3rk( M, N, NRHS, IOFFSET, NB, ABSTOL,
393 $ RELTOL, KP1, MAXC2NRM, A, LDA, DONE, KB,
394 $ MAXC2NRMK, RELMAXC2NRMK, JPIV, TAU,
395 $ VN1, VN2, AUXV, F, LDF, IWORK, INFO )
396 IMPLICIT NONE
397*
398* -- LAPACK auxiliary routine --
399* -- LAPACK is a software package provided by Univ. of Tennessee, --
400* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
401*
402* .. Scalar Arguments ..
403 LOGICAL DONE
404 INTEGER INFO, IOFFSET, KB, KP1, LDA, LDF, M, N,
405 $ NB, NRHS
406 REAL ABSTOL, MAXC2NRM, MAXC2NRMK, RELMAXC2NRMK,
407 $ reltol
408* ..
409* .. Array Arguments ..
410 INTEGER IWORK( * ), JPIV( * )
411 REAL VN1( * ), VN2( * )
412 COMPLEX A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * )
413* ..
414*
415* =====================================================================
416*
417* .. Parameters ..
418 REAL ZERO, ONE
419 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
420 COMPLEX CZERO, CONE
421 parameter( czero = ( 0.0e+0, 0.0e+0 ),
422 $ cone = ( 1.0e+0, 0.0e+0 ) )
423* ..
424* .. Local Scalars ..
425 INTEGER ITEMP, J, K, MINMNFACT, MINMNUPDT,
426 $ LSTICC, KP, I, IF
427 REAL HUGEVAL, TAUNAN, TEMP, TEMP2, TOL3Z
428 COMPLEX AIK
429* ..
430* .. External Subroutines ..
431 EXTERNAL cgemm, cgemv, clarfg, cswap
432* ..
433* .. Intrinsic Functions ..
434 INTRINSIC abs, real, conjg, imag, max, min, sqrt
435* ..
436* .. External Functions ..
437 LOGICAL SISNAN
438 INTEGER ISAMAX
439 REAL SLAMCH, SCNRM2
440 EXTERNAL sisnan, slamch, isamax, scnrm2
441* ..
442* .. Executable Statements ..
443*
444* Initialize INFO
445*
446 info = 0
447*
448* MINMNFACT in the smallest dimension of the submatrix
449* A(IOFFSET+1:M,1:N) to be factorized.
450*
451 minmnfact = min( m-ioffset, n )
452 minmnupdt = min( m-ioffset, n+nrhs )
453 nb = min( nb, minmnfact )
454 tol3z = sqrt( slamch( 'Epsilon' ) )
455 hugeval = slamch( 'Overflow' )
456*
457* Compute factorization in a while loop over NB columns,
458* K is the column index in the block A(1:M,1:N).
459*
460 k = 0
461 lsticc = 0
462 done = .false.
463*
464 DO WHILE ( k.LT.nb .AND. lsticc.EQ.0 )
465 k = k + 1
466 i = ioffset + k
467*
468 IF( i.EQ.1 ) THEN
469*
470* We are at the first column of the original whole matrix A_orig,
471* therefore we use the computed KP1 and MAXC2NRM from the
472* main routine.
473*
474 kp = kp1
475*
476 ELSE
477*
478* Determine the pivot column in K-th step, i.e. the index
479* of the column with the maximum 2-norm in the
480* submatrix A(I:M,K:N).
481*
482 kp = ( k-1 ) + isamax( n-k+1, vn1( k ), 1 )
483*
484* Determine the maximum column 2-norm and the relative maximum
485* column 2-norm of the submatrix A(I:M,K:N) in step K.
486*
487 maxc2nrmk = vn1( kp )
488*
489* ============================================================
490*
491* Check if the submatrix A(I:M,K:N) contains NaN, set
492* INFO parameter to the column number, where the first NaN
493* is found and return from the routine.
494* We need to check the condition only if the
495* column index (same as row index) of the original whole
496* matrix is larger than 1, since the condition for whole
497* original matrix is checked in the main routine.
498*
499 IF( sisnan( maxc2nrmk ) ) THEN
500*
501 done = .true.
502*
503* Set KB, the number of factorized partial columns
504* that are non-zero in each step in the block,
505* i.e. the rank of the factor R.
506* Set IF, the number of processed rows in the block, which
507* is the same as the number of processed rows in
508* the original whole matrix A_orig.
509*
510 kb = k - 1
511 IF = i - 1
512 info = kb + kp
513*
514* Set RELMAXC2NRMK to NaN.
515*
516 relmaxc2nrmk = maxc2nrmk
517*
518* There is no need to apply the block reflector to the
519* residual of the matrix A stored in A(KB+1:M,KB+1:N),
520* since the submatrix contains NaN and we stop
521* the computation.
522* But, we need to apply the block reflector to the residual
523* right hand sides stored in A(KB+1:M,N+1:N+NRHS), if the
524* residual right hand sides exist. This occurs
525* when ( NRHS != 0 AND KB <= (M-IOFFSET) ):
526*
527* A(I+1:M,N+1:N+NRHS) := A(I+1:M,N+1:N+NRHS) -
528* A(I+1:M,1:KB) * F(N+1:N+NRHS,1:KB)**H.
529
530 IF( nrhs.GT.0 .AND. kb.LT.(m-ioffset) ) THEN
531 CALL cgemm( 'No transpose', 'Conjugate transpose',
532 $ m-IF, nrhs, kb, -cone, a( if+1, 1 ), lda,
533 $ f( n+1, 1 ), ldf, cone, a( if+1, n+1 ), lda )
534 END IF
535*
536* There is no need to recompute the 2-norm of the
537* difficult columns, since we stop the factorization.
538*
539* Array TAU(KF+1:MINMNFACT) is not set and contains
540* undefined elements.
541*
542* Return from the routine.
543*
544 RETURN
545 END IF
546*
547* Quick return, if the submatrix A(I:M,K:N) is
548* a zero matrix. We need to check it only if the column index
549* (same as row index) is larger than 1, since the condition
550* for the whole original matrix A_orig is checked in the main
551* routine.
552*
553 IF( maxc2nrmk.EQ.zero ) THEN
554*
555 done = .true.
556*
557* Set KB, the number of factorized partial columns
558* that are non-zero in each step in the block,
559* i.e. the rank of the factor R.
560* Set IF, the number of processed rows in the block, which
561* is the same as the number of processed rows in
562* the original whole matrix A_orig.
563*
564 kb = k - 1
565 IF = i - 1
566 relmaxc2nrmk = zero
567*
568* There is no need to apply the block reflector to the
569* residual of the matrix A stored in A(KB+1:M,KB+1:N),
570* since the submatrix is zero and we stop the computation.
571* But, we need to apply the block reflector to the residual
572* right hand sides stored in A(KB+1:M,N+1:N+NRHS), if the
573* residual right hand sides exist. This occurs
574* when ( NRHS != 0 AND KB <= (M-IOFFSET) ):
575*
576* A(I+1:M,N+1:N+NRHS) := A(I+1:M,N+1:N+NRHS) -
577* A(I+1:M,1:KB) * F(N+1:N+NRHS,1:KB)**H.
578*
579 IF( nrhs.GT.0 .AND. kb.LT.(m-ioffset) ) THEN
580 CALL cgemm( 'No transpose', 'Conjugate transpose',
581 $ m-IF, nrhs, kb, -cone, a( if+1, 1 ), lda,
582 $ f( n+1, 1 ), ldf, cone, a( if+1, n+1 ), lda )
583 END IF
584*
585* There is no need to recompute the 2-norm of the
586* difficult columns, since we stop the factorization.
587*
588* Set TAUs corresponding to the columns that were not
589* factorized to ZERO, i.e. set TAU(KB+1:MINMNFACT) = CZERO,
590* which is equivalent to seting TAU(K:MINMNFACT) = CZERO.
591*
592 DO j = k, minmnfact
593 tau( j ) = czero
594 END DO
595*
596* Return from the routine.
597*
598 RETURN
599*
600 END IF
601*
602* ============================================================
603*
604* Check if the submatrix A(I:M,K:N) contains Inf,
605* set INFO parameter to the column number, where
606* the first Inf is found plus N, and continue
607* the computation.
608* We need to check the condition only if the
609* column index (same as row index) of the original whole
610* matrix is larger than 1, since the condition for whole
611* original matrix is checked in the main routine.
612*
613 IF( info.EQ.0 .AND. maxc2nrmk.GT.hugeval ) THEN
614 info = n + k - 1 + kp
615 END IF
616*
617* ============================================================
618*
619* Test for the second and third tolerance stopping criteria.
620* NOTE: There is no need to test for ABSTOL.GE.ZERO, since
621* MAXC2NRMK is non-negative. Similarly, there is no need
622* to test for RELTOL.GE.ZERO, since RELMAXC2NRMK is
623* non-negative.
624* We need to check the condition only if the
625* column index (same as row index) of the original whole
626* matrix is larger than 1, since the condition for whole
627* original matrix is checked in the main routine.
628*
629 relmaxc2nrmk = maxc2nrmk / maxc2nrm
630*
631 IF( maxc2nrmk.LE.abstol .OR. relmaxc2nrmk.LE.reltol ) THEN
632*
633 done = .true.
634*
635* Set KB, the number of factorized partial columns
636* that are non-zero in each step in the block,
637* i.e. the rank of the factor R.
638* Set IF, the number of processed rows in the block, which
639* is the same as the number of processed rows in
640* the original whole matrix A_orig;
641*
642 kb = k - 1
643 IF = i - 1
644*
645* Apply the block reflector to the residual of the
646* matrix A and the residual of the right hand sides B, if
647* the residual matrix and and/or the residual of the right
648* hand sides exist, i.e. if the submatrix
649* A(I+1:M,KB+1:N+NRHS) exists. This occurs when
650* KB < MINMNUPDT = min( M-IOFFSET, N+NRHS ):
651*
652* A(IF+1:M,K+1:N+NRHS) := A(IF+1:M,KB+1:N+NRHS) -
653* A(IF+1:M,1:KB) * F(KB+1:N+NRHS,1:KB)**H.
654*
655 IF( kb.LT.minmnupdt ) THEN
656 CALL cgemm( 'No transpose', 'Conjugate transpose',
657 $ m-IF, n+nrhs-kb, kb,-cone, a( if+1, 1 ), lda,
658 $ f( kb+1, 1 ), ldf, cone, a( if+1, kb+1 ), lda )
659 END IF
660*
661* There is no need to recompute the 2-norm of the
662* difficult columns, since we stop the factorization.
663*
664* Set TAUs corresponding to the columns that were not
665* factorized to ZERO, i.e. set TAU(KB+1:MINMNFACT) = CZERO,
666* which is equivalent to seting TAU(K:MINMNFACT) = CZERO.
667*
668 DO j = k, minmnfact
669 tau( j ) = czero
670 END DO
671*
672* Return from the routine.
673*
674 RETURN
675*
676 END IF
677*
678* ============================================================
679*
680* End ELSE of IF(I.EQ.1)
681*
682 END IF
683*
684* ===============================================================
685*
686* If the pivot column is not the first column of the
687* subblock A(1:M,K:N):
688* 1) swap the K-th column and the KP-th pivot column
689* in A(1:M,1:N);
690* 2) swap the K-th row and the KP-th row in F(1:N,1:K-1)
691* 3) copy the K-th element into the KP-th element of the partial
692* and exact 2-norm vectors VN1 and VN2. (Swap is not needed
693* for VN1 and VN2 since we use the element with the index
694* larger than K in the next loop step.)
695* 4) Save the pivot interchange with the indices relative to the
696* the original matrix A_orig, not the block A(1:M,1:N).
697*
698 IF( kp.NE.k ) THEN
699 CALL cswap( m, a( 1, kp ), 1, a( 1, k ), 1 )
700 CALL cswap( k-1, f( kp, 1 ), ldf, f( k, 1 ), ldf )
701 vn1( kp ) = vn1( k )
702 vn2( kp ) = vn2( k )
703 itemp = jpiv( kp )
704 jpiv( kp ) = jpiv( k )
705 jpiv( k ) = itemp
706 END IF
707*
708* Apply previous Householder reflectors to column K:
709* A(I:M,K) := A(I:M,K) - A(I:M,1:K-1)*F(K,1:K-1)**H.
710*
711 IF( k.GT.1 ) THEN
712 DO j = 1, k - 1
713 f( k, j ) = conjg( f( k, j ) )
714 END DO
715 CALL cgemv( 'No transpose', m-i+1, k-1, -cone, a( i, 1 ),
716 $ lda, f( k, 1 ), ldf, cone, a( i, k ), 1 )
717 DO j = 1, k - 1
718 f( k, j ) = conjg( f( k, j ) )
719 END DO
720 END IF
721*
722* Generate elementary reflector H(k) using the column A(I:M,K).
723*
724 IF( i.LT.m ) THEN
725 CALL clarfg( m-i+1, a( i, k ), a( i+1, k ), 1, tau( k ) )
726 ELSE
727 tau( k ) = czero
728 END IF
729*
730* Check if TAU(K) contains NaN, set INFO parameter
731* to the column number where NaN is found and return from
732* the routine.
733* NOTE: There is no need to check TAU(K) for Inf,
734* since CLARFG cannot produce TAU(KK) or Householder vector
735* below the diagonal containing Inf. Only BETA on the diagonal,
736* returned by CLARFG can contain Inf, which requires
737* TAU(K) to contain NaN. Therefore, this case of generating Inf
738* by CLARFG is covered by checking TAU(K) for NaN.
739*
740 IF( sisnan( real( tau(k) ) ) ) THEN
741 taunan = real( tau(k) )
742 ELSE IF( sisnan( imag( tau(k) ) ) ) THEN
743 taunan = imag( tau(k) )
744 ELSE
745 taunan = zero
746 END IF
747*
748 IF( sisnan( taunan ) ) THEN
749*
750 done = .true.
751*
752* Set KB, the number of factorized partial columns
753* that are non-zero in each step in the block,
754* i.e. the rank of the factor R.
755* Set IF, the number of processed rows in the block, which
756* is the same as the number of processed rows in
757* the original whole matrix A_orig.
758*
759 kb = k - 1
760 IF = i - 1
761 info = k
762*
763* Set MAXC2NRMK and RELMAXC2NRMK to NaN.
764*
765 maxc2nrmk = taunan
766 relmaxc2nrmk = taunan
767*
768* There is no need to apply the block reflector to the
769* residual of the matrix A stored in A(KB+1:M,KB+1:N),
770* since the submatrix contains NaN and we stop
771* the computation.
772* But, we need to apply the block reflector to the residual
773* right hand sides stored in A(KB+1:M,N+1:N+NRHS), if the
774* residual right hand sides exist. This occurs
775* when ( NRHS != 0 AND KB <= (M-IOFFSET) ):
776*
777* A(I+1:M,N+1:N+NRHS) := A(I+1:M,N+1:N+NRHS) -
778* A(I+1:M,1:KB) * F(N+1:N+NRHS,1:KB)**H.
779*
780 IF( nrhs.GT.0 .AND. kb.LT.(m-ioffset) ) THEN
781 CALL cgemm( 'No transpose', 'Conjugate transpose',
782 $ m-IF, nrhs, kb, -cone, a( if+1, 1 ), lda,
783 $ f( n+1, 1 ), ldf, cone, a( if+1, n+1 ), lda )
784 END IF
785*
786* There is no need to recompute the 2-norm of the
787* difficult columns, since we stop the factorization.
788*
789* Array TAU(KF+1:MINMNFACT) is not set and contains
790* undefined elements.
791*
792* Return from the routine.
793*
794 RETURN
795 END IF
796*
797* ===============================================================
798*
799 aik = a( i, k )
800 a( i, k ) = cone
801*
802* ===============================================================
803*
804* Compute the current K-th column of F:
805* 1) F(K+1:N,K) := tau(K) * A(I:M,K+1:N)**H * A(I:M,K).
806*
807 IF( k.LT.n+nrhs ) THEN
808 CALL cgemv( 'Conjugate transpose', m-i+1, n+nrhs-k,
809 $ tau( k ), a( i, k+1 ), lda, a( i, k ), 1,
810 $ czero, f( k+1, k ), 1 )
811 END IF
812*
813* 2) Zero out elements above and on the diagonal of the
814* column K in matrix F, i.e elements F(1:K,K).
815*
816 DO j = 1, k
817 f( j, k ) = czero
818 END DO
819*
820* 3) Incremental updating of the K-th column of F:
821* F(1:N,K) := F(1:N,K) - tau(K) * F(1:N,1:K-1) * A(I:M,1:K-1)**H
822* * A(I:M,K).
823*
824 IF( k.GT.1 ) THEN
825 CALL cgemv( 'Conjugate Transpose', m-i+1, k-1, -tau( k ),
826 $ a( i, 1 ), lda, a( i, k ), 1, czero,
827 $ auxv( 1 ), 1 )
828*
829 CALL cgemv( 'No transpose', n+nrhs, k-1, cone,
830 $ f( 1, 1 ), ldf, auxv( 1 ), 1, cone,
831 $ f( 1, k ), 1 )
832 END IF
833*
834* ===============================================================
835*
836* Update the current I-th row of A:
837* A(I,K+1:N+NRHS) := A(I,K+1:N+NRHS)
838* - A(I,1:K)*F(K+1:N+NRHS,1:K)**H.
839*
840 IF( k.LT.n+nrhs ) THEN
841 CALL cgemm( 'No transpose', 'Conjugate transpose',
842 $ 1, n+nrhs-k, k, -cone, a( i, 1 ), lda,
843 $ f( k+1, 1 ), ldf, cone, a( i, k+1 ), lda )
844 END IF
845*
846 a( i, k ) = aik
847*
848* Update the partial column 2-norms for the residual matrix,
849* only if the residual matrix A(I+1:M,K+1:N) exists, i.e.
850* when K < MINMNFACT = min( M-IOFFSET, N ).
851*
852 IF( k.LT.minmnfact ) THEN
853*
854 DO j = k + 1, n
855 IF( vn1( j ).NE.zero ) THEN
856*
857* NOTE: The following lines follow from the analysis in
858* Lapack Working Note 176.
859*
860 temp = abs( a( i, j ) ) / vn1( j )
861 temp = max( zero, ( one+temp )*( one-temp ) )
862 temp2 = temp*( vn1( j ) / vn2( j ) )**2
863 IF( temp2.LE.tol3z ) THEN
864*
865* At J-index, we have a difficult column for the
866* update of the 2-norm. Save the index of the previous
867* difficult column in IWORK(J-1).
868* NOTE: ILSTCC > 1, threfore we can use IWORK only
869* with N-1 elements, where the elements are
870* shifted by 1 to the left.
871*
872 iwork( j-1 ) = lsticc
873*
874* Set the index of the last difficult column LSTICC.
875*
876 lsticc = j
877*
878 ELSE
879 vn1( j ) = vn1( j )*sqrt( temp )
880 END IF
881 END IF
882 END DO
883*
884 END IF
885*
886* End of while loop.
887*
888 END DO
889*
890* Now, afler the loop:
891* Set KB, the number of factorized columns in the block;
892* Set IF, the number of processed rows in the block, which
893* is the same as the number of processed rows in
894* the original whole matrix A_orig, IF = IOFFSET + KB.
895*
896 kb = k
897 IF = i
898*
899* Apply the block reflector to the residual of the matrix A
900* and the residual of the right hand sides B, if the residual
901* matrix and and/or the residual of the right hand sides
902* exist, i.e. if the submatrix A(I+1:M,KB+1:N+NRHS) exists.
903* This occurs when KB < MINMNUPDT = min( M-IOFFSET, N+NRHS ):
904*
905* A(IF+1:M,K+1:N+NRHS) := A(IF+1:M,KB+1:N+NRHS) -
906* A(IF+1:M,1:KB) * F(KB+1:N+NRHS,1:KB)**H.
907*
908 IF( kb.LT.minmnupdt ) THEN
909 CALL cgemm( 'No transpose', 'Conjugate transpose',
910 $ m-IF, n+nrhs-kb, kb, -cone, a( if+1, 1 ), lda,
911 $ f( kb+1, 1 ), ldf, cone, a( if+1, kb+1 ), lda )
912 END IF
913*
914* Recompute the 2-norm of the difficult columns.
915* Loop over the index of the difficult columns from the largest
916* to the smallest index.
917*
918 DO WHILE( lsticc.GT.0 )
919*
920* LSTICC is the index of the last difficult column is greater
921* than 1.
922* ITEMP is the index of the previous difficult column.
923*
924 itemp = iwork( lsticc-1 )
925*
926* Compute the 2-norm explicilty for the last difficult column and
927* save it in the partial and exact 2-norm vectors VN1 and VN2.
928*
929* NOTE: The computation of VN1( LSTICC ) relies on the fact that
930* SCNRM2 does not fail on vectors with norm below the value of
931* SQRT(SLAMCH('S'))
932*
933 vn1( lsticc ) = scnrm2( m-IF, a( if+1, lsticc ), 1 )
934 vn2( lsticc ) = vn1( lsticc )
935*
936* Downdate the index of the last difficult column to
937* the index of the previous difficult column.
938*
939 lsticc = itemp
940*
941 END DO
942*
943 RETURN
944*
945* End of CLAQP3RK
946*
947 END
subroutine claqp3rk(m, n, nrhs, ioffset, nb, abstol, reltol, kp1, maxc2nrm, a, lda, done, kb, maxc2nrmk, relmaxc2nrmk, jpiv, tau, vn1, vn2, auxv, f, ldf, iwork, info)
CLAQP3RK computes a step of truncated QR factorization with column pivoting of a complex m-by-n matri...
Definition claqp3rk.f:396
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine clarfg(n, alpha, x, incx, tau)
CLARFG generates an elementary reflector (Householder matrix).
Definition clarfg.f:106
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81