Matrix Orthogonalization Routine

[yelu1220]

Hey guys,

The lab I work with is trying to port some of their old code that uses the NAG libraries over to a new supercomputer cluster that does not have NAG but has Intel MKL. Most of the general LAPACK and BLAS routines have a simple counterpart, but there is one NAG routine that I was wondering if MKL had something similar to.

The routine in NAG is called: F05AAF. It takes in a set of vectors and orthogonalizes them. We use it to orthogonalize a matrix in our code. The MKL libraries have orthogonal factorization routines, but I do not believe that it is the same thing as matrix orthogonalization (NAG uses the Gram-Schmidt process).

Here is the NAG description of F05AAF:http://www.nag.co.uk/numeric/Fl/manual/pdf/F05/f05aaf.pdf

Could anybody shine some light on whether MKL has this capability?

Thanks.

[Michael_C_Intel4]

Hello,

I'm an MKL LAPACK engineer. MKL doesn't have matrix orthogonalization routine. Generally, MKL supports LAPACK standard, which is described at www.netlib.org/lapack,but this routine isn't standard. Anyway, we can introduce functionality in MKL if it's useful for our customers.

Michael.

[Gennady_F_Intel]

Quoting - yelu1220

Hey guys,

The lab I work with is trying to port some of their old code that uses the NAG libraries over to a new supercomputer cluster that does not have NAG but has Intel MKL. Most of the general LAPACK and BLAS routines have a simple counterpart, but there is one NAG routine that I was wondering if MKL had something similar to.

The routine in NAG is called: F05AAF. It takes in a set of vectors and orthogonalizes them. We use it to orthogonalize a matrix in our code. The MKL libraries have orthogonal factorization routines, but I do not believe that it is the same thing as matrix orthogonalization (NAG uses the Gram-Schmidt process).

Here is the NAG description of F05AAF:http://www.nag.co.uk/numeric/Fl/manual/pdf/F05/f05aaf.pdf

Could anybody shine some light on whether MKL has this capability?

Thanks.

Yelu1220,

We would recommend you submit the issue against MKL to Premier support( https://premier.intel.com/ )

--Gennady

[Chao_Y_Intel]

Hi,

For such request, MKL LAPACK functionality can be used:

If we consider QR factorization of the M-by-N (M>N) matrx A=Q*R in the form where Q is square M-by-M matrix and R is an upper triangular M-by-N matrix. The equality A=Q*R can be re-written also as a product Q1*R1 where Q1 is a rectangular M-by-N submatrix of the matrix Q and R1 is M-by-M submatrix of the R. Let us note that columns of Q1 are orthonormal (they are orthogonal to each other and have norms equal to 1). The equality A=Q1*R1 can be treated as every column of A is a linear combination of Q1 columns, i.e. they span the same linear space. In other words, columns of Q1 is nothing else but a result of ortogonalization of columns A.

Of course, cases where R degenerates should be cosnsidered with some care and R must be trapezoidal if rank(A) is less than N.

QR functionality of MKL does not give Q directly- after calling DGEQRF a function DORGQR must be called.

DGEQPF might be useful for degenerate cases.

how does this works for you?

Thanks,
Chao

[Chao_Y_Intel]

Hi,

For such request, MKL LAPACK functionality can be used:

If we consider QR factorization of the M-by-N (M>N) matrx A=Q*R in the form where Q is square M-by-M matrix and R is an upper triangular M-by-N matrix. The equality A=Q*R can be re-written also as a product Q1*R1 where Q1 is a rectangular M-by-N submatrix of the matrix Q and R1 is M-by-M submatrix of the R. Let us note that columns of Q1 are orthonormal (they are orthogonal to each other and have norms equal to 1). The equality A=Q1*R1 can be treated as every column of A is a linear combination of Q1 columns, i.e. they span the same linear space. In other words, columns of Q1 is nothing else but a result of ortogonalization of columns A.

Of course, cases where R degenerates should be cosnsidered with some care and R must be trapezoidal if rank(A) is less than N.

QR functionality of MKL does not give Q directly- after calling DGEQRF a function DORGQR must be called.

DGEQPF might be useful for degenerate cases.

how does this works for you?

Thanks,
Chao