cpotrf.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 211 行 · 第 1/2 页
HTML
211 行
IF( .NOT.UPPER .AND. .NOT.<a name="LSAME.88"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -4
END IF
IF( INFO.NE.0 ) THEN
CALL <a name="XERBLA.96"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CPOTRF.96"></a><a href="cpotrf.f.html#CPOTRF.1">CPOTRF</a>'</span>, -INFO )
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span> IF( N.EQ.0 )
$ RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Determine the block size for this environment.
</span><span class="comment">*</span><span class="comment">
</span> NB = <a name="ILAENV.107"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="CPOTRF.107"></a><a href="cpotrf.f.html#CPOTRF.1">CPOTRF</a>'</span>, UPLO, N, -1, -1, -1 )
IF( NB.LE.1 .OR. NB.GE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Use unblocked code.
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="CPOTF2.112"></a><a href="cpotf2.f.html#CPOTF2.1">CPOTF2</a>( UPLO, N, A, LDA, INFO )
ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Use blocked code.
</span><span class="comment">*</span><span class="comment">
</span> IF( UPPER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute the Cholesky factorization A = U'*U.
</span><span class="comment">*</span><span class="comment">
</span> DO 10 J = 1, N, NB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Update and factorize the current diagonal block and test
</span><span class="comment">*</span><span class="comment"> for non-positive-definiteness.
</span><span class="comment">*</span><span class="comment">
</span> JB = MIN( NB, N-J+1 )
CALL CHERK( <span class="string">'Upper'</span>, <span class="string">'Conjugate transpose'</span>, JB, J-1,
$ -ONE, A( 1, J ), LDA, ONE, A( J, J ), LDA )
CALL <a name="CPOTF2.129"></a><a href="cpotf2.f.html#CPOTF2.1">CPOTF2</a>( <span class="string">'Upper'</span>, JB, A( J, J ), LDA, INFO )
IF( INFO.NE.0 )
$ GO TO 30
IF( J+JB.LE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute the current block row.
</span><span class="comment">*</span><span class="comment">
</span> CALL CGEMM( <span class="string">'Conjugate transpose'</span>, <span class="string">'No transpose'</span>, JB,
$ N-J-JB+1, J-1, -CONE, A( 1, J ), LDA,
$ A( 1, J+JB ), LDA, CONE, A( J, J+JB ),
$ LDA )
CALL CTRSM( <span class="string">'Left'</span>, <span class="string">'Upper'</span>, <span class="string">'Conjugate transpose'</span>,
$ <span class="string">'Non-unit'</span>, JB, N-J-JB+1, CONE, A( J, J ),
$ LDA, A( J, J+JB ), LDA )
END IF
10 CONTINUE
<span class="comment">*</span><span class="comment">
</span> ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute the Cholesky factorization A = L*L'.
</span><span class="comment">*</span><span class="comment">
</span> DO 20 J = 1, N, NB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Update and factorize the current diagonal block and test
</span><span class="comment">*</span><span class="comment"> for non-positive-definiteness.
</span><span class="comment">*</span><span class="comment">
</span> JB = MIN( NB, N-J+1 )
CALL CHERK( <span class="string">'Lower'</span>, <span class="string">'No transpose'</span>, JB, J-1, -ONE,
$ A( J, 1 ), LDA, ONE, A( J, J ), LDA )
CALL <a name="CPOTF2.158"></a><a href="cpotf2.f.html#CPOTF2.1">CPOTF2</a>( <span class="string">'Lower'</span>, JB, A( J, J ), LDA, INFO )
IF( INFO.NE.0 )
$ GO TO 30
IF( J+JB.LE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute the current block column.
</span><span class="comment">*</span><span class="comment">
</span> CALL CGEMM( <span class="string">'No transpose'</span>, <span class="string">'Conjugate transpose'</span>,
$ N-J-JB+1, JB, J-1, -CONE, A( J+JB, 1 ),
$ LDA, A( J, 1 ), LDA, CONE, A( J+JB, J ),
$ LDA )
CALL CTRSM( <span class="string">'Right'</span>, <span class="string">'Lower'</span>, <span class="string">'Conjugate transpose'</span>,
$ <span class="string">'Non-unit'</span>, N-J-JB+1, JB, CONE, A( J, J ),
$ LDA, A( J+JB, J ), LDA )
END IF
20 CONTINUE
END IF
END IF
GO TO 40
<span class="comment">*</span><span class="comment">
</span> 30 CONTINUE
INFO = INFO + J - 1
<span class="comment">*</span><span class="comment">
</span> 40 CONTINUE
RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> End of <a name="CPOTRF.184"></a><a href="cpotrf.f.html#CPOTRF.1">CPOTRF</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?