SUBROUTINECTRTRIF(UPLO,DIAG,N,A,LDA,INFO)**--LAPACKroutine(version3.0)--*Univ.ofTennessee,Univ.ofCaliforniaBerkeley,NAGLtd.,*CourantInstitute,ArgonneNationalLab,andRiceUniversity*September30,1994**..ScalarArguments..CHARACTERDIAG,UPLOINTEGERINFO,LDA,N*..*..ArrayArguments..COMPLEXA(LDA,*)*..**Purpose*=======**CTRTRIcomputestheinverseofacomplexupperorlowertriangular*matrixA.**ThisistheLevel3BLASversionofthealgorithm.**Arguments*=========**UPLO(input)CHARACTER*1*='U':Aisuppertriangular;*='L':Aislowertriangular.**DIAG(input)CHARACTER*1*='N':Aisnon-unittriangular;*='U':Aisunittriangular.**N(input)INTEGER*TheorderofthematrixA.N>=0.**A(input/output)COMPLEXarray,dimension(LDA,N)*Onentry,thetriangularmatrixA.IfUPLO='U',the*leadingN-by-NuppertriangularpartofthearrayAcontains*theuppertriangularmatrix,andthestrictlylower*triangularpartofAisnotreferenced.IfUPLO='L',the*leadingN-by-NlowertriangularpartofthearrayAcontains*thelowertriangularmatrix,andthestrictlyupper*triangularpartofAisnotreferenced.IfDIAG='U',the*diagonalelementsofAarealsonotreferencedandare*assumedtobe1.*Onexit,the(triangular)inverseoftheoriginalmatrix,in*thesamestorageformat.**LDA(input)INTEGER*TheleadingdimensionofthearrayA.LDA>=max(1,N).**INFO(output)INTEGER*=0:successfulexit*<0:ifINFO=-i,thei-thargumenthadanillegalvalue*>0:ifINFO=i,A(i,i)isexactlyzero.Thetriangular*matrixissingularanditsinversecannotbecomputed.**=====================================================================**..Parameters..COMPLEXONE,ZEROPARAMETER(ONE=(1.0E+0,0.0E+0),$ZERO=(0.0E+0,0.0E+0))*..*..LocalScalars..LOGICALNOUNIT,UPPERINTEGERJ,JB,NB,NN*..*..ExternalFunctions..LOGICALLSAMEEXTERNALLSAME*..*..ExternalSubroutines..EXTERNALCTRMM,CTRSM,CTRTI2,XERBLA*..*..IntrinsicFunctions..INTRINSICMAX,MIN*..*..ExecutableStatements..**Testtheinputparameters.*INFO=0UPPER=LSAME(UPLO,'U')NOUNIT=LSAME(DIAG,'N')IF(.NOT.UPPER.AND..NOT.LSAME(UPLO,'L'))THENINFO=-1ELSEIF(.NOT.NOUNIT.AND..NOT.LSAME(DIAG,'U'))THENINFO=-2ELSEIF(N.LT.0)THENINFO=-3ELSEIF(LDA.LT.MAX(1,N))THENINFO=-5ENDIFIF(INFO.NE.0)THENCALLXERBLA('CTRTRI',-INFO)RETURNENDIF**Quickreturnifpossible*IF(N.EQ.0)$RETURN**Checkforsingularityifnon-unit.*IF(NOUNIT)THENDO10INFO=1,NIF(A(INFO,INFO).EQ.ZERO)$RETURN10CONTINUEINFO=0ENDIF**Determinetheblocksizeforthisenvironment.*NB=128IF(NB.LE.1.OR.NB.GE.N)THEN**Useunblockedcode*CALLCTRTI2(UPLO,DIAG,N,A,LDA,INFO)ELSE**Useblockedcode*IF(UPPER)THEN**Computeinverseofuppertriangularmatrix*DO20J=1,N,NBJB=MIN(NB,N-J+1)**Computerows1:j-1ofcurrentblockcolumn*CALLCTRMM('Left','Upper','No transpose',DIAG,J-1,$JB,ONE,A,LDA,A(1,J),LDA)CALLCTRSM('Right','Upper','No transpose',DIAG,J-1,$JB,-ONE,A(J,J),LDA,A(1,J),LDA)**Computeinverseofcurrentdiagonalblock*CALLCTRTI2('Upper',DIAG,JB,A(J,J),LDA,INFO)20CONTINUEELSE**Computeinverseoflowertriangularmatrix*NN=((N-1)/NB)*NB+1DO30J=NN,1,-NBJB=MIN(NB,N-J+1)IF(J+JB.LE.N)THEN**Computerowsj+jb:nofcurrentblockcolumn*CALLCTRMM('Left','Lower','No transpose',DIAG,$N-J-JB+1,JB,ONE,A(J+JB,J+JB),LDA,$A(J+JB,J),LDA)CALLCTRSM('Right','Lower','No transpose',DIAG,$N-J-JB+1,JB,-ONE,A(J,J),LDA,$A(J+JB,J),LDA)ENDIF**Computeinverseofcurrentdiagonalblock*CALLCTRTI2('Lower',DIAG,JB,A(J,J),LDA,INFO)30CONTINUEENDIFENDIF*RETURN**EndofCTRTRI*END
此处可能存在不合适展示的内容,页面不予展示。您可通过相关编辑功能自查并修改。
如您确认内容无涉及 不当用语 / 纯广告导流 / 暴力 / 低俗色情 / 侵权 / 盗版 / 虚假 / 无价值内容或违法国家有关法律法规的内容,可点击提交进行申诉,我们将尽快为您处理。