Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
This community is designed for sharing of public information. Please do not share Intel or third-party confidential information here.

Reversed solve after zgetrf ?

Pierre_B_
Beginner
280 Views

Hello,

I'm trying to implement a blocked LDLT factorization but i'm facing some problems with the blocked column update.

Is there any way to solve the equation X.A=B instead of A.X=B (zgetrs) after a LU decomposition of A with zgetrf ?

Thank you in advance,

Pierre

0 Kudos
1 Solution
Konstantin_A_Intel
280 Views

Hi Pierre,

That's possible. First, after ZGETRF we have A=P*L*U (where P is row permutation storing in IPIV array).

So, X*A=B  <=>  X*P*L*U=B  <=>  ((X*P)*L)*U=B.  Let Y=X*P and  Z=Y*L.  In this case, you can solve:

(1) Z*U = B  - find Z using ZTRSM routine

(2) Y*L = Z - find Y using ZTRSM routine

(3) X*P = Y  <=> X=Y*P^-1  -  find X applying inverse permutation P.

Regards,

Konstantin

View solution in original post

3 Replies
Konstantin_A_Intel
281 Views

Hi Pierre,

That's possible. First, after ZGETRF we have A=P*L*U (where P is row permutation storing in IPIV array).

So, X*A=B  <=>  X*P*L*U=B  <=>  ((X*P)*L)*U=B.  Let Y=X*P and  Z=Y*L.  In this case, you can solve:

(1) Z*U = B  - find Z using ZTRSM routine

(2) Y*L = Z - find Y using ZTRSM routine

(3) X*P = Y  <=> X=Y*P^-1  -  find X applying inverse permutation P.

Regards,

Konstantin

Pierre_B_
Beginner
280 Views

Hello Konstantin,

Thank you very much for your answer. I did not know i could use ztrsm on the output of zgetrf.

What i am supposed to do to apply the inverse permutation ? Is there a routine for that or what should I do ?

Thank you,

Pierre

mecej4
Black Belt
280 Views

The simplest solution, if you can rearrange your work so that you factorize AT instead of A, is to note that X A = B is the same as ATXT = BT. The last equation can be solved by a single call to ZGESV.

If, however, you wish to use the steps of #2, ...:

To undo the permutations, you can do the corresponding row interchanges after obtaining Y in order to obtain X, or you can form and use the inverse of the permutation vector when printing the solution. Here is an example (needs more testing for correctness!). For input data, you can use the file sgetrsx.d in the MKL examples/lapack/data directory.

! Example program to illustrate solving X.A = B
!    1. Call ?GETRF to factorize A = P L U
!    2. Call ?TRSM  to solve Z U = B
!    3. Call ?TRSM  to solve Y L = Z
!    4. Apply inverse of P to retrieve X = Y inv(P)
!
! Ref.: software.intel.com/en-us/forums/intel-math-kernel-library/topic/702864
!
    program strsmx
       implicit none
       integer nin, nout
       parameter (nin=5, nout=6)
       integer nmax, lda, nrhmax, ldb
       parameter (nmax=8, lda=nmax, nrhmax=nmax, ldb=nrhmax)
!     .. Local Variables ..
       integer i, ifail, info, j, k, n, nrhs
       real a(lda, nmax), b(nrhmax, nmax), alpha
       integer ipiv(nmax), iipiv(nmax)
!     .. Executable Statements ..
       write (nout, *) 'STRSM Example Program Results'
!     Skip heading in data file
       read (nin, *)
       read (nin, *) n, nrhs
       if (n<=nmax .and. nrhs<=nrhmax) then
!
!        Read A and B from data file
!
          read (nin, *)((a(i,j),j=1,n), i=1, n)
          read (nin, *)((b(j,i),j=1,nrhs), i=1, n)
!
!        Factorize A
!
          call sgetrf(n, n, a, lda, ipiv, info)
!
          write (nout, *)
          if (info/=0) then
             write (nout, *) 'The matrix A is singular'
          end if
!
!           Compute solution
!
          alpha = 1.0
          call strsm('R', 'U', 'N', 'N', n, n, alpha, a, lda, b, ldb)
          call strsm('R', 'L', 'N', 'U', n, n, alpha, a, lda, b, ldb)
!
!        Form inverse of pivot array
!
          iipiv(1:n) = (/ (i,i=1,n) /)
          do i = n, 1, -1
             j = ipiv(i)
             if (j/=i) then
                k = iipiv(i)
                iipiv(i) = iipiv(j)
                iipiv(j) = k
             end if
          end do
!
!           Print solution
!
          do i = 1, n
             write (*, '(10(2x,ES12.3))')(b(j,iipiv(i)), j=1, nrhs)
          end do
       end if
       stop
!
    end program strsmx

 

Reply