1，源码
SUBROUTINE INIT_MATRIX(A, m, n, lda)
DOUBLE PRECISION A(*)
CALL SRAND(2024)
DO i=1, m
DO j=1, n
A(i + lda*(j-1)) = RAND() + RAND()
C WRITE(*, '(F8.4)') A(i)
END DO
END DO
END
SUBROUTINE PRINT_MATRIX(A, m, n, lda)
DOUBLE PRECISION A(*)
DO i=1, m
DO j=1, n
C WRITE(*, '(F8.4 , $)') A(i + lda*(j-1))
WRITE(*, 20) A(i + lda*(j-1))
20 FORMAT (1X,F8.4, $)
END DO
PRINT *,''
END DO
END
PROGRAM MATRIX_EX
INTEGER m, n, lda
DOUBLE PRECISION A(20)
lda = 4
m = 4
n = 5
CALL INIT_MATRIX(A, m, n, lda)
PRINT *, "A ="
CALL PRINT_MATRIX(A, m, n, lda)
PRINT *, "A(11) ="
PRINT *, A(11)
ENDMakefile
matrix_ex: matrix_ex.f
gfortran -g $< -o $@
.PHONY: clean
clean:
-rm matrix_ex2，执行
make
./matrix_ex效果:
3, dgesvd example
* 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.
* =============================================================================
*
* DGESVD Example.
* ==============
*
* Program computes the singular value decomposition of a general
* rectangular matrix A:
*
* 8.79 9.93 9.83 5.45 3.16
* 6.11 6.91 5.04 -0.27 7.98
* -9.15 -7.93 4.86 4.85 3.01
* 9.57 1.64 8.83 0.74 5.80
* -3.49 4.02 9.80 10.00 4.27
* 9.84 0.15 -8.99 -6.02 -5.31
*
* Description.
* ============
*
* The routine computes the singular value decomposition (SVD) of a real
* m-by-n matrix A, optionally computing the left and/or right singular
* vectors. The SVD is written as
*
* A = U*SIGMA*VT
*
* 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 orthogonal matrix and VT (V transposed)
* is an n-by-n orthogonal 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 VT, not V.
*
* Example Program Results.
* ========================
*
* DGESVD Example Program Results
*
* Singular values
* 27.47 22.64 8.56 5.99 2.01
*
* Left singular vectors (stored columnwise)
* -0.59 0.26 0.36 0.31 0.23
* -0.40 0.24 -0.22 -0.75 -0.36
* -0.03 -0.60 -0.45 0.23 -0.31
* -0.43 0.24 -0.69 0.33 0.16
* -0.47 -0.35 0.39 0.16 -0.52
* 0.29 0.58 -0.02 0.38 -0.65
*
* Right singular vectors (stored rowwise)
* -0.25 -0.40 -0.69 -0.37 -0.41
* 0.81 0.36 -0.25 -0.37 -0.10
* -0.26 0.70 -0.22 0.39 -0.49
* 0.40 -0.45 0.25 0.43 -0.62
* -0.22 0.14 0.59 -0.63 -0.44
* =============================================================================
*
* .. Parameters ..
INTEGER M, N
PARAMETER ( M = 6, N = 5 )
INTEGER LDA, LDU, LDVT
PARAMETER ( LDA = M, LDU = M, LDVT = N )
INTEGER LWMAX
PARAMETER ( LWMAX = 1000 )
*
* .. Local Scalars ..
INTEGER INFO, LWORK
*
* .. Local Arrays ..
DOUBLE PRECISION A( LDA, N ), U( LDU, M ), VT( LDVT, N ), S( N ),
$ WORK( LWMAX )
DATA A/
$ 8.79, 6.11,-9.15, 9.57,-3.49, 9.84,
$ 9.93, 6.91,-7.93, 1.64, 4.02, 0.15,
$ 9.83, 5.04, 4.86, 8.83, 9.80,-8.99,
$ 5.45,-0.27, 4.85, 0.74,10.00,-6.02,
$ 3.16, 7.98, 3.01, 5.80, 4.27,-5.31
$ /
*
* .. External Subroutines ..
EXTERNAL DGESVD
EXTERNAL PRINT_MATRIX
*
* .. Intrinsic Functions ..
INTRINSIC INT, MIN
*
* .. Executable Statements ..
WRITE(*,*)'DGESVD Example Program Results'
*
* Query the optimal workspace.
*
LWORK = -1
CALL DGESVD( 'All', 'All', M, N, A, LDA, S, U, LDU, VT, LDVT,
$ WORK, LWORK, INFO )
LWORK = MIN( LWMAX, INT( WORK( 1 ) ) )
*
* Compute SVD.
*
CALL DGESVD( 'All', 'All', M, N, A, LDA, S, U, LDU, VT, LDVT,
$ WORK, LWORK, 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_MATRIX( 'Singular values', 1, N, S, 1 )
*
* Print left singular vectors.
*
CALL PRINT_MATRIX( 'Left singular vectors (stored columnwise)',
$ M, N, U, LDU )
*
* Print right singular vectors.
*
CALL PRINT_MATRIX( 'Right singular vectors (stored rowwise)',
$ N, N, VT, LDVT )
STOP
END
*
* End of DGESVD Example.
*
* =============================================================================
*
* Auxiliary routine: printing a matrix.
*
SUBROUTINE PRINT_MATRIX( 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
Makefile
svd_dgesve_ex: svd_dgesve_ex.f
gfortran -g $< ../lapack-3.11/liblapack.a ../lapack-3.11/librefblas.a -o $@运行
