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