実一般矩形行列のQR分解

C言語によるサンプルソースコード : 使用関数名:nag_dgeqrf (f08aec)

ホーム > 製品 > nAG数値計算ライブラリ > サンプルソースコード集 > 実一般矩形行列のQR分解 (C言語/C++)

Keyword: 実一般矩形行列, QR分解

概要

本サンプルは実一般矩形行列のQR分解を行うC言語によるサンプルプログラムです。 本サンプルは以下に示される最小二乗問題をQR分解を行って解き、解を出力します。

実一般矩形行のデータ 

※本サンプルはnAG Cライブラリに含まれる関数 nag_dgeqrf() のExampleコードです。本サンプル及び関数の詳細情報は nag_dgeqrf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

(本関数の詳細はnag_dgeqrf のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11
12
13
14

このデータをダウンロード
nag_dgeqrf (f08aec) 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
 -3.15   2.19
 -0.11  -3.64
  1.99   0.57
 -2.70   8.23
  0.26  -6.35
  4.50  -1.48                 :End of matrix B

  • 1行目はタイトル行で読み飛ばされます。
  • 2行目に行列Aの行数(m)、列数(n)、右辺の数(nrhs)を指定しています。
  • 3~8行目に行列Aの要素を指定しています。
  • 9~14行目に行列Bの要素を指定しています。

出力結果

(本関数の詳細はnag_dgeqrf のマニュアルページを参照)
1
2
3
4
5
6
7
8

この出力例をダウンロード
nag_dgeqrf (f08aec) Example Program Results

 Least squares solution(s)
            1          2
 1      1.5146    -1.5838
 2      1.8621     0.5536
 3     -1.4467     1.3491
 4      0.0396     2.9600

  • 5~9行目に最小二乗解が出力されています。

ソースコード

(本関数の詳細はnag_dgeqrf のマニュアルページを参照)

※本サンプルソースコードはnAG数値計算ライブラリ(Windows, Linux, MAC等に対応)の関数を呼び出します。
サンプルのコンパイル及び実行方法

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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

このソースコードをダウンロード
/* nag_dgeqrf (f08aec) Example Program.
 *
 * CLL6I261D/CLL6I261DL Version.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagf07.h>
#include <nagf08.h>
#include <nagx04.h>

int main(void)
{
  /* Scalars */
  Integer i, j, m, n, nrhs, pda, pdb, tau_len;
  Integer exit_status = 0;
  NagError fail;
  Nag_OrderType order;
  /* Arrays */
  double *a = 0, *b = 0, *tau = 0;
#ifdef nAG_LOAD_FP
  /* The following line is needed to force the Microsoft linker
     to load floating point support */
  float force_loading_of_ms_float_support = 0;
#endif /* nAG_LOAD_FP */

#ifdef nAG_COLUMN_MAJOR
#define A(I, J) a[(J - 1) * pda + I - 1]
#define B(I, J) b[(J - 1) * pdb + I - 1]
  order = Nag_ColMajor;
#else
#define A(I, J) a[(I - 1) * pda + J - 1]
#define B(I, J) b[(I - 1) * pdb + J - 1]
  order = Nag_RowMajor;
#endif

  INIT_FAIL(fail);

  printf("nag_dgeqrf (f08aec) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  scanf("%ld%ld%ld%*[^\n] ", &m, &n, &nrhs);
#ifdef nAG_COLUMN_MAJOR
  pda = m;
  pdb = m;
#else
  pda = n;
  pdb = nrhs;
#endif
  tau_len = MIN(m, n);

  /* Allocate memory */
  if (!(a = nAG_ALLOC(m * n, double)) ||
      !(b = nAG_ALLOC(m * nrhs, double)) ||
      !(tau = nAG_ALLOC(tau_len, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  /* Read A and B from data file */
  for (i = 1; i <= m; ++i) {
    for (j = 1; j <= n; ++j)
      scanf("%lf", &A(i, j));
  }
  scanf("%*[^\n] ");
  for (i = 1; i <= m; ++i) {
    for (j = 1; j <= nrhs; ++j)
      scanf("%lf", &B(i, j));
  }
  scanf("%*[^\n] ");

  /* Compute the QR factorization of A */
  /* nag_dgeqrf (f08aec).
   * QR factorization of real general rectangular matrix
   */
  nag_dgeqrf(order, m, n, a, pda, tau, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_dgeqrf (f08aec).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Compute C = (Q^T)*B, storing the result in B */
  /* nag_dormqr (f08agc).
   * Apply orthogonal transformation determined by nag_dgeqrf (f08aec)
   * or nag_dgeqpf (f08bec)
   */
  nag_dormqr(order, Nag_LeftSide, Nag_Trans, m, nrhs, n, a, pda,
             tau, b, pdb, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_dormqr (f08agc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Compute least squares solution by back-substitution in R*X = C */
  /* nag_dtrtrs (f07tec).
   * Solution of real triangular system of linear equations,
   * multiple right-hand sides
   */
  nag_dtrtrs(order, Nag_Upper, Nag_NoTrans, Nag_NonUnitDiag, n, nrhs,
             a, pda, b, pdb, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_dtrtrs (f07tec).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  /* Print least squares solution(s) */
  /* nag_gen_real_mat_print (x04cac).
   * Print real general matrix (easy-to-use)
   */
  fflush(stdout);
  nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, nrhs,
                         b, pdb, "Least squares solution(s)", 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
END:
  nAG_FREE(a);
  nAG_FREE(b);
  nAG_FREE(tau);
  return exit_status;
}


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