Visible to Intel only — GUID: GUID-63EE5070-EBA2-4D7A-AD3A-22C0388825DA
ZGESDD Example Program in Fortran
* Copyright (C) 2009-2015 Intel Corporation. All Rights Reserved.
* The information and material ("Material") provided below is owned by Intel
* Corporation or its suppliers or licensors, and title to such Material remains
* with Intel Corporation or its suppliers or licensors. The Material contains
* proprietary information of Intel or its suppliers and licensors. The Material
* is protected by worldwide copyright laws and treaty provisions. No part of
* the Material may be copied, reproduced, published, uploaded, posted,
* transmitted, or distributed in any way without Intel's prior express written
* permission. No license under any patent, copyright or other intellectual
* property rights in the Material is granted to or conferred upon you, either
* expressly, by implication, inducement, estoppel or otherwise. Any license
* under such intellectual property rights must be express and approved by Intel
* in writing.
* =============================================================================
*
* ZGESDD Example.
* ==============
*
* Program computes the singular value decomposition of a general
* rectangular complex matrix A using a divide and conquer method, where A is:
*
* ( -5.40, 7.40) ( 6.00, 6.38) ( 9.91, 0.16) ( -5.28, -4.16)
* ( 1.09, 1.55) ( 2.60, 0.07) ( 3.98, -5.26) ( 2.03, 1.11)
* ( 9.88, 1.91) ( 4.92, 6.31) ( -2.11, 7.39) ( -9.81, -8.98)
*
* Description.
* ============
*
* The routine computes the singular value decomposition (SVD) of a complex
* m-by-n matrix A, optionally computing the left and/or right singular
* vectors. If singular vectors are desired, it uses a divide and conquer
* algorithm. The SVD is written as
*
* A = U*SIGMA*VH
*
* where SIGMA is an m-by-n matrix which is zero except for its min(m,n)
* diagonal elements, U is an m-by-m unitary matrix and VH (V conjugate
* transposed) is an n-by-n unitary matrix. The diagonal elements of SIGMA
* are the singular values of A; they are real and non-negative, and are
* returned in descending order. The first min(m, n) columns of U and V are
* the left and right singular vectors of A.
*
* Note that the routine returns VH, not V.
*
* Example Program Results.
* ========================
*
* ZGESDD Example Program Results
*
* Singular values
* 21.76 16.60 3.97
*
* Left singular vectors (stored columnwise)
* ( 0.55, 0.00) ( 0.76, 0.00) ( -0.34, 0.00)
* ( -0.04, -0.15) ( 0.27, -0.23) ( 0.55, -0.74)
* ( 0.81, 0.12) ( -0.52, -0.14) ( 0.13, -0.11)
*
* Right singular vectors (stored rowwise)
* ( 0.23, 0.21) ( 0.37, 0.39) ( 0.24, 0.33) ( -0.56, -0.37)
* ( -0.58, 0.40) ( 0.11, 0.17) ( 0.60, -0.27) ( 0.16, 0.06)
* ( 0.60, 0.12) ( -0.19, 0.30) ( 0.39, 0.20) ( 0.45, 0.31)
* =============================================================================
*
* .. Parameters ..
INTEGER M, N
PARAMETER ( M = 3, N = 4 )
INTEGER LDA, LDU, LDVT
PARAMETER ( LDA = M, LDU = M, LDVT = N )
INTEGER LWMAX
PARAMETER ( LWMAX = 1000 )
*
* .. Local Scalars ..
INTEGER INFO, LWORK
*
* .. Local Arrays ..
* IWORK dimension should be at least 8*MIN(M,N)
INTEGER IWORK( 8*M )
* RWORK dimension should be at least 5*(MIN(M,N))**2 + 7*MIN(M,N))
DOUBLE PRECISION S( M ), RWORK( 5*M*M + 7*M )
COMPLEX*16 A( LDA, N ), U( LDU, M ), VT( LDVT, N ),
$ WORK( LWMAX )
DATA A/
$ (-5.40, 7.40),( 1.09, 1.55),( 9.88, 1.91),
$ ( 6.00, 6.38),( 2.60, 0.07),( 4.92, 6.31),
$ ( 9.91, 0.16),( 3.98,-5.26),(-2.11, 7.39),
$ (-5.28,-4.16),( 2.03, 1.11),(-9.81,-8.98)
$ /
*
* .. External Subroutines ..
EXTERNAL ZGESDD
EXTERNAL PRINT_MATRIX, PRINT_RMATRIX
*
* .. Intrinsic Functions ..
INTRINSIC INT, MIN
*
* .. Executable Statements ..
WRITE(*,*)'ZGESDD Example Program Results'
*
* Query the optimal workspace.
*
LWORK = -1
CALL ZGESDD( 'Singular vectors', M, N, A, LDA, S, U, LDU, VT,
$ LDVT, WORK, LWORK, RWORK, IWORK, INFO )
LWORK = MIN( LWMAX, INT( WORK( 1 ) ) )
*
* Compute SVD.
*
CALL ZGESDD( 'Singular vectors', M, N, A, LDA, S, U, LDU, VT,
$ LDVT, WORK, LWORK, RWORK, IWORK, INFO )
*
* Check for convergence.
*
IF( INFO.GT.0 ) THEN
WRITE(*,*)'The algorithm computing SVD failed to converge.'
STOP
END IF
*
* Print singular values.
*
CALL PRINT_RMATRIX( 'Singular values', 1, M, S, 1 )
*
* Print left singular vectors.
*
CALL PRINT_MATRIX( 'Left singular vectors (stored columnwise)',
$ M, M, U, LDU )
*
* Print right singular vectors.
*
CALL PRINT_MATRIX( 'Right singular vectors (stored rowwise)',
$ M, N, VT, LDVT )
STOP
END
*
* End of ZGESDD Example.
*
* =============================================================================
*
* Auxiliary routine: printing a matrix.
*
SUBROUTINE PRINT_MATRIX( DESC, M, N, A, LDA )
CHARACTER*(*) DESC
INTEGER M, N, LDA
COMPLEX*16 A( LDA, * )
*
INTEGER I, J
*
WRITE(*,*)
WRITE(*,*) DESC
DO I = 1, M
WRITE(*,9998) ( A( I, J ), J = 1, N )
END DO
*
9998 FORMAT( 11(:,1X,'(',F6.2,',',F6.2,')') )
RETURN
END
*
* Auxiliary routine: printing a real matrix.
*
SUBROUTINE PRINT_RMATRIX( DESC, M, N, A, LDA )
CHARACTER*(*) DESC
INTEGER M, N, LDA
DOUBLE PRECISION A( LDA, * )
*
INTEGER I, J
*
WRITE(*,*)
WRITE(*,*) DESC
DO I = 1, M
WRITE(*,9998) ( A( I, J ), J = 1, N )
END DO
*
9998 FORMAT( 11(:,1X,F6.2) )
RETURN
END
Parent topic: ZGESDD Example