Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- Reversed solve after zgetrf ?

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Pierre_B_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-09-2016
07:41 AM

131 Views

Reversed solve after zgetrf ?

Hello,

I'm trying to implement a blocked LDL^{T} 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

1 Solution

Konstantin_A_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-10-2016
08:43 PM

131 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

Link Copied

3 Replies

Konstantin_A_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-10-2016
08:43 PM

132 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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-11-2016
08:59 AM

131 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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-12-2016
05:57 AM

131 Views

The simplest solution, if you can rearrange your work so that you factorize A^{T} instead of A, is to note that X A = B is the same as A^{T}X^{T} = B^{T}. 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

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.