ORTH
Orthogonalizes the columns of a matrix.
Function Return Value
Orthogonal matrix Q from the decomposition A=QR. If A is rank-3, Q is rank-3. (Output)
Required Argument
A — Matrix A to be decomposed. Must be an array of rank-2 or rank-3 (box data) of type real, double, complex, or double complex. (Input)
Optional Arguments, Packaged Options
R — Upper-triangular or upper trapezoidal matrix R, from the QR decomposition. If this optional argument is present, the decomposition is complete. If A is rank-3, R is rank-3. (Output)
The option and derived type names are given in the following tables:
|
Option Name for ORTH |
Option Value |
|
Skip_error_processing |
5 |
|
Name of Unallocated Option Array to Use for Setting Options |
Use |
Derived Type |
|
?_orth_options(:) |
Use when setting options for calls hereafter. |
?_options |
|
?_orth_options_once(:) |
Use when setting options for next call only. |
?_options |
For a description on how to use these options, see Matrix Optional Data Changes.
FORTRAN 90 Interface
ORTH (A [, …])
Description
Orthogonalizes the columns of a matrix. The decomposition A = QR is computed using a forward and backward sweep of the Modified Gram-Schmidt algorithm.
Examples
Example 1: (Operator_ex19.f90)
use linear_operators
use lin_sol_tri_int
use rand_int
use Numerical_Libraries
implicit none
! This is the equivalent of Example 3 (using operators) for LIN_SOL_TRI.
integer i, nopt
integer, parameter :: n=128, k=n/4, ncoda=1, lda=2
real(kind(1e0)), parameter :: s_one=1e0, s_zero=0e0
real(kind(1e0)) A(lda,n), EVAL(k)
type(s_options) :: iopt(2)
real(kind(1e0)) d(n), b(n), d_t(2*n,k), c_t(2*n,k), perf_ratio, &
b_t(2*n,k), y_t(2*n,k), eval_t(k), res(n,k)
logical small
! This flag is used to get the k largest eigenvalues.
small = .false.
! Generate the main diagonal and the co-diagonal of the
! tridiagonal matrix.
b=rand(b); d=rand(d)
A(1,1:)=b; A(2,1:)=d
! Use Numerical Libraries routine for the calculation of k
! largest eigenvalues.
CALL EVASB (N, K, A, LDA, NCODA, SMALL, EVAL)
EVAL_T = EVAL
! Use IMSL Fortran Numerical Librarytridiagonal solver for inverse iteration
! calculation of eigenvectors.
factorization_choice: do nopt=0,1
! Create k tridiagonal problems, one for each inverse
! iteration system.
b_t(1:n,1:k) = spread(b,DIM=2,NCOPIES=k)
c_t(1:n,1:k) = EOSHIFT(b_t(1:n,1:k),SHIFT=1,DIM=1)
d_t(1:n,1:k) = spread(d,DIM=2,NCOPIES=k) - &
spread(EVAL_T,DIM=1,NCOPIES=n)
! Start the right-hand side at random values, scaled downward
! to account for the expected 'blowup' in the solution.
y_t=rand(y_t)
! Do two iterations for the eigenvectors.
do i=1, 2
y_t(1:n,1:k) = y_t(1:n,1:k)*epsilon(s_one)
call lin_sol_tri(c_t, d_t, b_t, y_t, &
iopt=iopt)
iopt(nopt+1) = s_lin_sol_tri_solve_only
end do
! Orthogonalize the eigenvectors. (This is the most
! intensive part of the computing.)
y_t(1:n,1:k) = ORTH(y_t(1:n,1:k))
! See if the performance ratio is smaller than the value one.
! If it is not the code will re-solve the systems using Gaussian
! Elimination. This is an exceptional event. It is a necessary
! complication for achieving reliable results.
res(1:n,1:k) = spread(d,DIM=2,NCOPIES=k)*y_t(1:n,1:k) + &
spread(b,DIM=2,NCOPIES=k)* &
EOSHIFT(y_t(1:n,1:k),SHIFT=-1,DIM=1) + &
EOSHIFT(spread(b,DIM=2,NCOPIES=k)*y_t(1:n,1:k),SHIFT=1) &
- y_t(1:n,1:k)*spread(EVAL_T(1:k),DIM=1,NCOPIES=n)
! If the factorization method is Cyclic Reduction and perf_ratio is
! larger than one, re-solve using Gaussian Elimination. If the
! method is already Gaussian Elimination, the loop exits
! and perf_ratio is checked at the end.
perf_ratio = norm(res(1:n,1:k),1) / &
norm(EVAL_T(1:k),1) / &
epsilon(s_one) / (5*n)
if (perf_ratio <= s_one) exit factorization_choice
iopt(nopt+1) = s_lin_sol_tri_use_Gauss_elim
end do factorization_choice
if (perf_ratio <= s_one) then
write (*,*) 'Example 3 for LIN_SOL_TRI (operators) is correct.'
end if
end
use linear_operators
use mpi_setup_int
integer, parameter :: N=32, nr=4
real (kind(1.e0)) A(N,N,nr), Q(N,N,nr)
! Setup for MPI
mp_nprocs = mp_setup()
if (mp_rank == 0) then
A = rand(A)
end if
Q = orth(A)
mp_nprocs = mp_setup ('Final')
end
