dlalsa.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 387 行 · 第 1/2 页
HTML
387 行
</span> INFO = 0
<span class="comment">*</span><span class="comment">
</span> IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
INFO = -1
ELSE IF( SMLSIZ.LT.3 ) THEN
INFO = -2
ELSE IF( N.LT.SMLSIZ ) THEN
INFO = -3
ELSE IF( NRHS.LT.1 ) THEN
INFO = -4
ELSE IF( LDB.LT.N ) THEN
INFO = -6
ELSE IF( LDBX.LT.N ) THEN
INFO = -8
ELSE IF( LDU.LT.N ) THEN
INFO = -10
ELSE IF( LDGCOL.LT.N ) THEN
INFO = -19
END IF
IF( INFO.NE.0 ) THEN
CALL <a name="XERBLA.196"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLALSA.196"></a><a href="dlalsa.f.html#DLALSA.1">DLALSA</a>'</span>, -INFO )
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Book-keeping and setting up the computation tree.
</span><span class="comment">*</span><span class="comment">
</span> INODE = 1
NDIML = INODE + N
NDIMR = NDIML + N
<span class="comment">*</span><span class="comment">
</span> CALL <a name="DLASDT.206"></a><a href="dlasdt.f.html#DLASDT.1">DLASDT</a>( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
$ IWORK( NDIMR ), SMLSIZ )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The following code applies back the left singular vector factors.
</span><span class="comment">*</span><span class="comment"> For applying back the right singular vector factors, go to 50.
</span><span class="comment">*</span><span class="comment">
</span> IF( ICOMPQ.EQ.1 ) THEN
GO TO 50
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The nodes on the bottom level of the tree were solved
</span><span class="comment">*</span><span class="comment"> by <a name="DLASDQ.217"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>. The corresponding left and right singular vector
</span><span class="comment">*</span><span class="comment"> matrices are in explicit form. First apply back the left
</span><span class="comment">*</span><span class="comment"> singular vector matrices.
</span><span class="comment">*</span><span class="comment">
</span> NDB1 = ( ND+1 ) / 2
DO 10 I = NDB1, ND
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IC : center row of each node
</span><span class="comment">*</span><span class="comment"> NL : number of rows of left subproblem
</span><span class="comment">*</span><span class="comment"> NR : number of rows of right subproblem
</span><span class="comment">*</span><span class="comment"> NLF: starting row of the left subproblem
</span><span class="comment">*</span><span class="comment"> NRF: starting row of the right subproblem
</span><span class="comment">*</span><span class="comment">
</span> I1 = I - 1
IC = IWORK( INODE+I1 )
NL = IWORK( NDIML+I1 )
NR = IWORK( NDIMR+I1 )
NLF = IC - NL
NRF = IC + 1
CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
$ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
$ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
10 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Next copy the rows of B that correspond to unchanged rows
</span><span class="comment">*</span><span class="comment"> in the bidiagonal matrix to BX.
</span><span class="comment">*</span><span class="comment">
</span> DO 20 I = 1, ND
IC = IWORK( INODE+I-1 )
CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
20 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Finally go through the left singular vector matrices of all
</span><span class="comment">*</span><span class="comment"> the other subproblems bottom-up on the tree.
</span><span class="comment">*</span><span class="comment">
</span> J = 2**NLVL
SQRE = 0
<span class="comment">*</span><span class="comment">
</span> DO 40 LVL = NLVL, 1, -1
LVL2 = 2*LVL - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> find the first node LF and last node LL on
</span><span class="comment">*</span><span class="comment"> the current level LVL
</span><span class="comment">*</span><span class="comment">
</span> IF( LVL.EQ.1 ) THEN
LF = 1
LL = 1
ELSE
LF = 2**( LVL-1 )
LL = 2*LF - 1
END IF
DO 30 I = LF, LL
IM1 = I - 1
IC = IWORK( INODE+IM1 )
NL = IWORK( NDIML+IM1 )
NR = IWORK( NDIMR+IM1 )
NLF = IC - NL
NRF = IC + 1
J = J - 1
CALL <a name="DLALS0.277"></a><a href="dlals0.f.html#DLALS0.1">DLALS0</a>( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
$ B( NLF, 1 ), LDB, PERM( NLF, LVL ),
$ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
$ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
$ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
$ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
$ INFO )
30 CONTINUE
40 CONTINUE
GO TO 90
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ICOMPQ = 1: applying back the right singular vector factors.
</span><span class="comment">*</span><span class="comment">
</span> 50 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> First now go through the right singular vector matrices of all
</span><span class="comment">*</span><span class="comment"> the tree nodes top-down.
</span><span class="comment">*</span><span class="comment">
</span> J = 0
DO 70 LVL = 1, NLVL
LVL2 = 2*LVL - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Find the first node LF and last node LL on
</span><span class="comment">*</span><span class="comment"> the current level LVL.
</span><span class="comment">*</span><span class="comment">
</span> IF( LVL.EQ.1 ) THEN
LF = 1
LL = 1
ELSE
LF = 2**( LVL-1 )
LL = 2*LF - 1
END IF
DO 60 I = LL, LF, -1
IM1 = I - 1
IC = IWORK( INODE+IM1 )
NL = IWORK( NDIML+IM1 )
NR = IWORK( NDIMR+IM1 )
NLF = IC - NL
NRF = IC + 1
IF( I.EQ.LL ) THEN
SQRE = 0
ELSE
SQRE = 1
END IF
J = J + 1
CALL <a name="DLALS0.322"></a><a href="dlals0.f.html#DLALS0.1">DLALS0</a>( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
$ BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
$ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
$ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
$ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
$ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
$ INFO )
60 CONTINUE
70 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The nodes on the bottom level of the tree were solved
</span><span class="comment">*</span><span class="comment"> by <a name="DLASDQ.333"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>. The corresponding right singular vector
</span><span class="comment">*</span><span class="comment"> matrices are in explicit form. Apply them back.
</span><span class="comment">*</span><span class="comment">
</span> NDB1 = ( ND+1 ) / 2
DO 80 I = NDB1, ND
I1 = I - 1
IC = IWORK( INODE+I1 )
NL = IWORK( NDIML+I1 )
NR = IWORK( NDIMR+I1 )
NLP1 = NL + 1
IF( I.EQ.ND ) THEN
NRP1 = NR
ELSE
NRP1 = NR + 1
END IF
NLF = IC - NL
NRF = IC + 1
CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
$ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
$ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
80 CONTINUE
<span class="comment">*</span><span class="comment">
</span> 90 CONTINUE
<span class="comment">*</span><span class="comment">
</span> RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> End of <a name="DLALSA.360"></a><a href="dlalsa.f.html#DLALSA.1">DLALSA</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?