複素一般化線形最小二乗: GLM

LAPACKサンプルソースコード : 使用ルーチン名:ZGGGLM

概要

本サンプルはFortran言語によりLAPACKルーチンZGGGLMを利用するサンプルプログラムです。

加重最小二乗問題を解きます。

\begin{displaymath}
\min_{x} \left\Vert B^{-1} (d - A x) \right\Vert _{2}
\end{displaymath}

ここで

\begin{displaymath}
A = \left(
\begin{array}{rrr}
0.96 - 0.81 i & -0.03 + 0.9...
...28 i & 0.20 - 0.12 i & -0.07 + 1.23 i
\end{array} \right), \:
\end{displaymath}


\begin{displaymath}
B = \left(
\begin{array}{cccc}
0.5 - 1.0 i & 0 & 0 & 0 \\...
... - 3.0 i & 0 \\
0 & 0 & 0 & 5.0 - 4.0 i
\end{array} \right)
\end{displaymath}

及び

\begin{displaymath}
d = \left(
\begin{array}{r}
6.00 - 0.40 i \\
-5.27 + 0.90 i \\
2.72 - 2.13 i \\
-1.30 - 2.80 i
\end{array} \right).
\end{displaymath}

入力データ

(本ルーチンの詳細はZGGGLM のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

このデータをダウンロード
ZGGGLM Example Program Data

  4             3             4                         :Values of M, N and P

( 0.96,-0.81) (-0.03, 0.96) (-0.91, 2.06)
(-0.98, 1.98) (-1.20, 0.19) (-0.66, 0.42)
( 0.62,-0.46) ( 1.01, 0.02) ( 0.63,-0.17)
( 1.08,-0.28) ( 0.20,-0.12) (-0.07, 1.23)               :End of matrix A

( 0.50,-1.00) ( 0.00, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)
( 0.00, 0.00) ( 1.00,-2.00) ( 0.00, 0.00) ( 0.00, 0.00)
( 0.00, 0.00) ( 0.00, 0.00) ( 2.00,-3.00) ( 0.00, 0.00)
( 0.00, 0.00) ( 0.00, 0.00) ( 0.00, 0.00) ( 5.00,-4.00) :End of matrix B

( 6.00,-0.40)
(-5.27, 0.90)
( 2.72,-2.13)
(-1.30,-2.80)                                           :End of vector d

出力結果

(本ルーチンの詳細はZGGGLM のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11

この出力例をダウンロード
 ZGGGLM Example Program Results

 Weighted least squares solution
 (  -0.9846,   1.9950) (   3.9929,  -4.9748) (  -3.0026,   0.9994)

 Residual vector
 ( 1.26E-04,-4.66E-04) ( 1.11E-03,-8.61E-04) ( 3.84E-03,-1.82E-03)
 ( 2.03E-03, 3.02E-03)

 Square root of the residual sum of squares
   5.79E-03

ソースコード

(本ルーチンの詳細はZGGGLM のマニュアルページを参照)

※本サンプルソースコードのご利用手順は「サンプルのコンパイル及び実行方法」をご参照下さい。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

このソースコードをダウンロード
    Program zggglm_example

!     ZGGGLM Example Program Text

!     Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com

!     .. Use Statements ..
      Use blas_interfaces, Only: dznrm2
      Use lapack_interfaces, Only: zggglm
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: rnorm
      Integer :: i, info, lda, ldb, lwork, m, n, p
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), b(:, :), d(:), work(:), x(:), &
        y(:)
!     .. Executable Statements ..
      Write (nout, *) 'ZGGGLM Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n, p
      lda = m
      ldb = m
      lwork = n + m + nb*(m+p)
      Allocate (a(lda,n), b(ldb,p), d(m), work(lwork), x(n), y(p))

!     Read A, B and D from data file

      Read (nin, *)(a(i,1:n), i=1, m)
      Read (nin, *)(b(i,1:p), i=1, m)
      Read (nin, *) d(1:m)

!     Solve the weighted least squares problem

!     minimize ||inv(B)*(d - A*x)|| (in the 2-norm)

      Call zggglm(m, n, p, a, lda, b, ldb, d, x, y, work, lwork, info)

!     Print least squares solution

      Write (nout, *) 'Weighted least squares solution'
      Write (nout, 100) x(1:n)

!     Print residual vector y = inv(B)*(d - A*x)

      Write (nout, *)
      Write (nout, *) 'Residual vector'
      Write (nout, 110) y(1:p)

!     Compute and print the square root of the residual sum of squares
      rnorm = dznrm2(p, y, 1)

      Write (nout, *)
      Write (nout, *) 'Square root of the residual sum of squares'
      Write (nout, 120) rnorm

100   Format (3(' (',F9.4,',',F9.4,')',:))
110   Format (3(' (',1P,E9.2,',',1P,E9.2,')',:))
120   Format (1X, 1P, E10.2)
    End Program


ご案内
関連情報
© 日本ニューメリカルアルゴリズムズグループ株式会社 2025
Privacy Policy  /  Trademarks