Program f01gbfe ! F01GBF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: dgemm, f01gbf, nag_wp, x04caf ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: t, tr Integer :: i, ifail, irevcm, lda, ldb, ldb2, & ldx, ldy, m, n ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: a(:,:), b(:,:), b2(:,:), comm(:), & p(:), r(:), x(:,:), y(:,:), z(:) Integer, Allocatable :: icomm(:) ! .. Executable Statements .. Write (nout,*) 'F01GBF 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 memory Allocate (a(lda,n)) Allocate (b(ldb,m)) Allocate (b2(ldb2,m)) Allocate (comm(n*m+3*n+12)) Allocate (x(ldx,2)) Allocate (y(ldy,2)) Allocate (icomm(2*n+40)) 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 Do i = 1, n tr = tr + a(i,i) End Do ! Find exp(tA)B irevcm = 0 ifail = 0 ! Initial call to f01gbf reverse communication interface revcm: Do Call f01gbf(irevcm,n,m,b,ldb,t,tr,b2,ldb2,x,ldx,y,ldy,p,r,z,comm, & icomm,ifail) If (irevcm==0) Then Exit revcm Else If (irevcm==1) Then ! Compute AB and store in B2 Call dgemm('N','N',n,m,n,1.0_nag_wp,a,lda,b,ldb,0.0_nag_wp,b2,ldb2) Else If (irevcm==2) Then ! Compute AX and store in Y Call dgemm('N','N',n,2,n,1.0_nag_wp,a,lda,x,ldx,0.0_nag_wp,y,ldy) Else If (irevcm==3) Then ! Compute A^T Y and store in X Call dgemm('T','N',n,2,n,1.0_nag_wp,a,lda,y,ldy,0.0_nag_wp,x,ldx) Else If (irevcm==4) Then ! Compute AZ and store in P Call dgemm('N','N',n,1,n,1.0_nag_wp,a,lda,z,n,0.0_nag_wp,p,n) Else ! Compute A^T Z and store in R Call dgemm('T','N',n,1,n,1.0_nag_wp,a,lda,z,n,0.0_nag_wp,r,n) End If End Do revcm If (ifail==0) Then ! Print solution ifail = 0 Call x04caf('G','N',n,m,b,ldb,'exp(tA)B',ifail) End If End Program f01gbfe