Program f01hbfe ! F01HBF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: f01hbf, nag_wp, x04daf, zgemm ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Complex (Kind=nag_wp) :: t, tr Integer :: i, ifail, irevcm, lda, ldb, ldb2, & ldx, ldy, m, n ! .. Local Arrays .. Complex (Kind=nag_wp), Allocatable :: a(:,:), b(:,:), b2(:,:), ccomm(:), & p(:), r(:), x(:,:), y(:,:), z(:) Real (Kind=nag_wp), Allocatable :: comm(:) Integer, Allocatable :: icomm(:) ! .. Executable Statements .. Write (nout,*) 'F01HBF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) Read (nin,*) n, m, t lda = n ldb = n ldb2 = n ldx = n ldy = n ! Allocate required space Allocate (a(lda,n)) Allocate (b(ldb,m)) Allocate (b2(ldb2,m)) Allocate (ccomm(n*(m+2))) Allocate (x(ldx,2)) Allocate (y(ldy,2)) Allocate (icomm(2*n+40)) Allocate (comm(14+3*n)) Allocate (p(n)) Allocate (r(n)) Allocate (z(n)) ! Read A from data file Read (nin,*)(a(i,1:n),i=1,n) ! Read B from data file Read (nin,*)(b(i,1:m),i=1,n) ! Compute the trace of A tr = (0.0_nag_wp,0.0_nag_wp) Do i = 1, n tr = tr + a(i,i) End Do ! Find exp(tA)B irevcm = 0 ifail = 0 ! Initial call to reverse communication interface f01hbf revcm: Do Call f01hbf(irevcm,n,m,b,ldb,t,tr,b2,ldb2,x,ldx,y,ldy,p,r,z,ccomm, & comm,icomm,ifail) If (irevcm==0) Then Exit revcm Else If (irevcm==1) Then ! Compute AB and store in B2 Call zgemm('N','N',n,m,n,(1.0_nag_wp,0.0_nag_wp),a,lda,b,ldb, & (0.0_nag_wp,0.0_nag_wp),b2,ldb2) Else If (irevcm==2) Then ! Compute AX and store in Y Call zgemm('N','N',n,2,n,(1.0_nag_wp,0.0_nag_wp),a,lda,x,ldx, & (0.0_nag_wp,0.0_nag_wp),y,ldy) Else If (irevcm==3) Then ! Compute A^H Y and store in X Call zgemm('C','N',n,2,n,(1.0_nag_wp,0.0_nag_wp),a,lda,y,ldy, & (0.0_nag_wp,0.0_nag_wp),x,ldx) Else If (irevcm==4) Then ! Compute Az and store in p Call zgemm('N','N',n,1,n,(1.0_nag_wp,0.0_nag_wp),a,lda,z,n, & (0.0_nag_wp,0.0_nag_wp),p,n) Else If (irevcm==5) Then ! Compute A^H z and store in r Call zgemm('C','N',n,1,n,(1.0_nag_wp,0.0_nag_wp),a,lda,z,n, & (0.0_nag_wp,0.0_nag_wp),r,n) End If ! Return to f01hbf End Do revcm If (ifail==0) Then ! Print solution ifail = 0 Call x04daf('G','N',n,m,b,ldb,'exp(tA)B',ifail) End If End Program f01hbfe