Keyword: 実一般矩形行列, QR分解
概要
本サンプルは実一般矩形行列のQR分解を行うFortranによるサンプルプログラムです。 本サンプルは以下に示される最小二乗問題をQR分解を行って解き、解を出力します。
※本サンプルはnAG Fortranライブラリに含まれるルーチン f08aef() のExampleコードです。本サンプル及びルーチンの詳細情報は f08aef のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はf08aef のマニュアルページを参照）| このデータをダウンロード |
F08AEF Example Program Data 6 4 2 :Values of M, N and NRHS -0.57 -1.28 -0.39 0.25 -1.93 1.08 -0.31 -2.14 2.30 0.24 0.40 -0.35 -1.93 0.64 -0.66 0.08 0.15 0.30 0.15 -2.13 -0.02 1.03 -1.43 0.50 :End of matrix A -2.67 0.41 -0.55 -3.10 3.34 -4.01 -0.77 2.76 0.48 -6.17 4.10 0.21 :End of matrix B
- １行目はタイトル行で読み飛ばされます。
- ３行目に行列Aの行数(m)、列数(n)、右辺の数(nrhs)を指定しています。
- ５〜１０行目に行列Aの要素を指定しています。
- １２〜１７行目に行列Bの要素を指定しています。
出力結果
（本ルーチンの詳細はf08aef のマニュアルページを参照）| この出力例をダウンロード |
F08AEF Example Program Results
Least-squares solution(s)
1 2
1 1.5339 -1.5753
2 1.8707 0.5559
3 -1.5241 1.3119
4 0.0392 2.9585
Square root(s) of the residual sum(s) of squares
2.22E-02 1.38E-02
- ５〜８行目に最小二乗解が出力されています。
- １１行目に残差平方和の平方根が出力されています。
ソースコード
（本ルーチンの詳細はf08aef のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
PROGRAM f08aefe
! F08AEF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
! .. Use Statements ..
USE nag_library, ONLY : dgeqrf, dnrm2, dormqr, dtrtrs, nag_wp, x04caf
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nb = 64, nin = 5, nout = 6
! .. Local Scalars ..
INTEGER :: i, ifail, info, j, lda, ldb, lwork, &
m, n, nrhs
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), b(:,:), rnorm(:), tau(:), &
work(:)
! .. Executable Statements ..
WRITE (nout,*) 'F08AEF Example Program Results'
WRITE (nout,*)
FLUSH (nout)
! Skip heading in data file
READ (nin,*)
READ (nin,*) m, n, nrhs
lda = m
ldb = m
lwork = nb*n
ALLOCATE (a(lda,n),b(ldb,nrhs),rnorm(nrhs),tau(n),work(lwork))
! Read A and B from data file
READ (nin,*) (a(i,1:n),i=1,m)
READ (nin,*) (b(i,1:nrhs),i=1,m)
! Compute the QR factorization of A
! The nAG name equivalent of dgeqrf is f08aef
CALL dgeqrf(m,n,a,lda,tau,work,lwork,info)
! Compute C = (C1) = (Q**T)*B, storing the result in B
! (C2)
! The nAG name equivalent of dormqr is f08agf
CALL dormqr('Left','Transpose',m,nrhs,n,a,lda,tau,b,ldb,work,lwork, &
info)
! Compute least-squares solutions by backsubstitution in
! R*X = C1
! The nAG name equivalent of dtrtrs is f07tef
CALL dtrtrs('Upper','No transpose','Non-Unit',n,nrhs,a,lda,b,ldb,info)
IF (info>0) THEN
WRITE (nout,*) 'The upper triangular factor, R, of A is singular, '
WRITE (nout,*) 'the least squares solution could not be computed'
ELSE
! Print least-squares solutions
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
CALL x04caf('General',' ',n,nrhs,b,ldb,'Least-squares solution(s)', &
ifail)
! Compute and print estimates of the square roots of the residual
! sums of squares
! The nAG name equivalent of dnrm2 is f06ejf
DO j = 1, nrhs
rnorm(j) = dnrm2(m-n,b(n+1,j),1)
END DO
WRITE (nout,*)
WRITE (nout,*) 'Square root(s) of the residual sum(s) of squares'
WRITE (nout,99999) rnorm(1:nrhs)
END IF
99999 FORMAT (5X,1P,7E11.2)
END PROGRAM f08aefe