Visible to Intel only — GUID: GUID-A601B7DC-A51A-48C5-B16B-5928F1ABFEEE
ZGELSD 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.
* =============================================================================
*
* ZGELSD Example.
* ==============
*
* Program computes the minimum norm-solution to a complex linear least squares
* problem using the singular value decomposition of A,
* where A is the coefficient matrix:
*
* ( 4.55, -0.32) ( -4.36, -4.76) ( 3.99, -6.84) ( 8.03, -6.47)
* ( 8.87, -3.11) ( 0.02, 8.43) ( 5.43, -9.30) ( 2.28, 8.94)
* ( -0.74, 1.16) ( 3.80, -6.12) ( -7.24, 0.72) ( 2.21, 9.52)
*
* and B is the right-hand side matrix:
*
* ( -8.25, 7.98) ( 2.91, -8.81)
* ( -5.04, 3.33) ( 6.19, 0.19)
* ( 7.98, -4.38) ( -5.96, 7.18)
*
* Description.
* ============
*
* The routine computes the minimum-norm solution to a complex linear least
* squares problem: minimize ||b - A*x|| using the singular value
* decomposition (SVD) of A. A is an m-by-n matrix which may be rank-deficient.
*
* Several right hand side vectors b and solution vectors x can be handled
* in a single call; they are stored as the columns of the m-by-nrhs right
* hand side matrix B and the n-by-nrhs solution matrix X.
*
* The effective rank of A is determined by treating as zero those singular
* values which are less than rcond times the largest singular value.
*
* Example Program Results.
* ========================
*
* ZGELSD Example Program Results
*
* Minimum norm solution
* ( -0.08, 0.09) ( 0.04, 0.16)
* ( -0.17, 0.10) ( 0.17, -0.47)
* ( -0.92, -0.01) ( 0.71, -0.41)
* ( -0.47, -0.26) ( 0.69, 0.02)
*
* Effective rank = 3
*
* Singular values
* 20.01 18.21 7.88
* =============================================================================
*
* .. Parameters ..
INTEGER M, N, NRHS
PARAMETER ( M = 3, N = 4, NRHS = 2 )
INTEGER LDA, LDB
PARAMETER ( LDA = M, LDB = N )
INTEGER LWMAX
PARAMETER ( LWMAX = 1000 )
*
* .. Local Scalars ..
INTEGER INFO, LWORK, RANK
DOUBLE PRECISION RCOND
*
* .. Local Arrays ..
* IWORK dimension should be at least 3*MIN(M,N)*NLVL + 11*MIN(M,N),
* RWORK dimension should be at least 10*MIN(M,N)+2*MIN(M,N)*SMLSIZ+
* +8*MIN(M,N)*NLVL+3*SMLSIZ*NRHS+(SMLSIZ+1)**2,
* where NLVL = MAX( 0, INT( LOG_2( MIN(M,N)/(SMLSIZ+1) ) )+1 )
* and SMLSIZ = 25
INTEGER IWORK( 3*M*0+11*M )
DOUBLE PRECISION S( M ), RWORK( 10*M+2*M*25+8*M*0+3*25*NRHS+26*26
$ )
COMPLEX*16 A( LDA, N ), B( LDB, NRHS ), WORK( LWMAX )
DATA A/
$ ( 4.55,-0.32),( 8.87,-3.11),(-0.74, 1.16),
$ (-4.36,-4.76),( 0.02, 8.43),( 3.80,-6.12),
$ ( 3.99,-6.84),( 5.43,-9.30),(-7.24, 0.72),
$ ( 8.03,-6.47),( 2.28, 8.94),( 2.21, 9.52)
$ /
DATA B/
$ (-8.25, 7.98),(-5.04, 3.33),( 7.98,-4.38),( 0.00, 0.00),
$ ( 2.91,-8.81),( 6.19, 0.19),(-5.96, 7.18),( 0.00, 0.00)
$ /
*
* .. External Subroutines ..
EXTERNAL ZGELSD
EXTERNAL PRINT_MATRIX, PRINT_RMATRIX
*
* .. Intrinsic Functions ..
INTRINSIC INT, MIN
*
* .. Executable Statements ..
WRITE(*,*)'ZGELSD Example Program Results'
* Negative RCOND means using default (machine precision) value
RCOND = -1.0
*
* Query the optimal workspace.
*
LWORK = -1
CALL ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK,
$ LWORK, RWORK, IWORK, INFO )
LWORK = MIN( LWMAX, INT( WORK( 1 ) ) )
*
* Solve the equations A*X = B.
*
CALL ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK,
$ LWORK, RWORK, IWORK, INFO )
*
* Check for convergence.
*
IF( INFO.GT.0 ) THEN
WRITE(*,*)'The algorithm computing SVD failed to converge;'
WRITE(*,*)'the least squares solution could not be computed.'
STOP
END IF
*
* Print minimum norm solution.
*
CALL PRINT_MATRIX( 'Minimum norm solution', N, NRHS, B, LDB )
*
* Print effective rank.
*
WRITE(*,'(/A,I6)')' Effective rank = ', RANK
*
* Print singular values.
*
CALL PRINT_RMATRIX( 'Singular values', 1, M, S, 1 )
STOP
END
*
* End of ZGELSD 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: ZGELSD Example