実行列ペアの一般化特異値分解

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

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

Keyword: 実行列, ペア, 一般化特異値分解

概要

本サンプルは実行列ペアの一般化特異値分解を行うC言語によるサンプルプログラムです。 本サンプルは以下に示される実行列ペアの一般化特異値分解を行い、条件数の推定値と一般化特異値の誤差推定値を出力します。

実行列のデータ 

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

入力データ

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

このデータをダウンロード
nag_dggsvd (f08vac) Example Program Data

  4    3    2   :Values of M, N and P

  1.0  2.0  3.0
  3.0  2.0  1.0
  4.0  5.0  6.0
  7.0  8.0  8.0 :End of matrix A

 -2.0 -3.0  3.0
  4.0  6.0  5.0 :End of matrix B

  • 1行目はタイトル行で読み飛ばされます。
  • 3行目に行列Aの行数(m)、行列Aと行列Bの列数(n)、行列Bの行数(p)を指定しています。
  • 5~8行目に行列Aの要素を指定しています。
  • 10~11行目に行列Bの要素を指定しています。

出力結果

(本関数の詳細はnag_dggsvd のマニュアルページを参照)
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

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

Number of infinite generalized singular values (k)
    1

Number of finite generalized singular values (l)
    2
Numerical rank of ( A^T B^T)^T   (k+l)
    3

Finite generalized singular values
   1.3151e+00    8.0185e-02

 Orthogonal matrix U
               1             2             3             4
 1   -1.3484e-01    5.2524e-01   -2.0924e-01    8.1373e-01
 2    6.7420e-01   -5.2213e-01   -3.8886e-01    3.4874e-01
 3    2.6968e-01    5.2757e-01   -6.5782e-01   -4.6499e-01
 4    6.7420e-01    4.1615e-01    6.1014e-01   -1.7776e-15

 Orthogonal matrix V
               1             2
 1    3.5539e-01   -9.3472e-01
 2    9.3472e-01    3.5539e-01

 Orthogonal matrix Q
               1             2             3
 1   -8.3205e-01   -9.4633e-02   -5.4657e-01
 2    5.5470e-01   -1.4195e-01   -8.1985e-01
 3    0.0000e+00   -9.8534e-01    1.7060e-01

 Nonsingular upper triangular matrix R
               1             2             3
 1   -2.0569e+00   -9.0121e+00   -9.3705e+00
 2                 -1.0882e+01   -7.2688e+00
 3                               -6.0405e+00

Estimate of reciprocal condition number for R
    4.2e-02

Error estimate for the generalized singular values
    2.6e-15

  • 4行目に無限大の一般化特異値の数が出力されています。
  • 7行目に有限の一般化特異値の数が出力されています。
  • 9行目に行列の有効ランクが出力されています。
  • 12行目に有限の一般化特異値が出力されています。
  • 16~20行目に直交行列Uが出力されています。
  • 24~26行目に直交行列Vが出力されています。
  • 30~33行目に直交行列Qが出力されています。
  • 37~40行目に正則上三角行列Rが出力されています。
  • 43行目にRの条件数の逆数の推定値が出力されています。
  • 46行目に一般化特異値の誤差推定が出力されています。

ソースコード

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

※本サンプルソースコードは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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215

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

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

int main(void)
{

  /* Scalars */
  double d, eps, rcond, serrbd;
  Integer exit_status = 0, i, irank, j, k, l, m, n, p,
         pda, pdb, pdq, pdu, pdv;
  NagError fail;
  Nag_OrderType order;

  /* Arrays */
  double *a = 0, *alpha = 0, *b = 0, *beta = 0, *q = 0, *u = 0, *v = 0;
  Integer *iwork = 0;

#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_dggsvd (f08vac) Example Program Results\n\n");
  /* Skip heading in data file */
  scanf("%*[^\n] ");

  scanf("%ld%ld%ld%*[^\n] ", &m, &n, &p);
  if (m <= 10 && n <= 10 && p <= 10) {
    /* Allocate memory */
    if (!(a = nAG_ALLOC(m * n, double)) ||
        !(alpha = nAG_ALLOC(n, double)) ||
        !(b = nAG_ALLOC(p * n, double)) ||
        !(beta = nAG_ALLOC(n, double)) ||
        !(q = nAG_ALLOC(n * n, double)) ||
        !(u = nAG_ALLOC(m * m, double)) ||
        !(v = nAG_ALLOC(p * p, double)) || !(iwork = nAG_ALLOC(n, Integer)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
#ifdef nAG_COLUMN_MAJOR
    pda = m;
    pdb = p;
    pdq = n;
    pdu = m;
    pdv = p;
#else
    pda = n;
    pdb = n;
    pdq = n;
    pdu = m;
    pdv = p;
#endif
  }
  else {
    printf("m and/or n too small\n");
    goto END;
  }
  /* Read the m by n matrix A and p by n matrix 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 <= p; ++i)
    for (j = 1; j <= n; ++j)
      scanf("%lf", &B(i, j));
  scanf("%*[^\n] ");

  /* nag_dggsvd (f08vac)
   * Compute the generalized singular value decomposition of (A, B)
   * (A = U*D1*(0 R)*(Q^T), B = V*D2*(0 R)*(Q^T), m.ge.n)
   */

  nag_dggsvd(order, Nag_AllU, Nag_ComputeV, Nag_ComputeQ, m, n, p, &k, &l, a,
             pda, b, pdb, alpha, beta, u, pdu, v, pdv, q, pdq, iwork, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_dggsvd (f08vac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  /* Print solution */
  irank = k + l;

  printf("Number of infinite generalized singular values (k)\n");
  printf("%5ld\n", k);
  printf("\nNumber of finite generalized singular values (l)\n");
  printf("%5ld\n", l);
  printf("Numerical rank of ( A^T B^T)^T   (k+l)\n");
  printf("%5ld\n", irank);

  printf("\nFinite generalized singular values\n");
  for (j = k; j < irank; ++j) {
    d = alpha[j] / beta[j];
    printf("%13.4e%s", d, (j + 1) % 8 == 0 || (j + 1) == irank ? "\n" : " ");
  }

  printf("\n");
  fflush(stdout);
  nag_gen_real_mat_print_comp(order, Nag_GeneralMatrix, Nag_NonUnitDiag, m,
                              m, u, pdu, "%13.4e", "Orthogonal matrix U",
                              Nag_IntegerLabels, 0, Nag_IntegerLabels, 0,
                              80, 0, 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\n");
  fflush(stdout);
  nag_gen_real_mat_print_comp(order, Nag_GeneralMatrix, Nag_NonUnitDiag, p,
                              p, v, pdv, "%13.4e", "Orthogonal matrix V",
                              Nag_IntegerLabels, 0, Nag_IntegerLabels, 0,
                              80, 0, 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\n");
  fflush(stdout);
  nag_gen_real_mat_print_comp(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n,
                              n, q, pdq, "%13.4e", "Orthogonal matrix Q",
                              Nag_IntegerLabels, 0, Nag_IntegerLabels, 0,
                              80, 0, 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\n");
  fflush(stdout);
  nag_gen_real_mat_print_comp(order, Nag_UpperMatrix, Nag_NonUnitDiag, irank,
                              irank, &A(1, n - irank + 1), pda, "%13.4e",
                              "Nonsingular upper triangular matrix R",
                              Nag_IntegerLabels, 0, Nag_IntegerLabels, 0,
                              80, 0, 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* nag_dtrcon (f07tgc)
   * estimate the reciprocal condition number of R
   */

  nag_dtrcon(order, Nag_InfNorm, Nag_Upper, Nag_NonUnitDiag, irank,
             &A(1, n - irank + 1), pda, &rcond, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_dtrcon (f07tgc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\nEstimate of reciprocal condition number for R\n");
  printf("%11.1e\n\n", rcond);

  /* So long as irank = n, get the machine precision, eps, and compute the
   * approximate error bound for the computed generalized singular values
   */
  if (irank == n) {
    eps = nag_machine_precision;
    serrbd = eps / rcond;
    printf("Error estimate for the generalized singular values");
    printf("\n%11.1e\n", serrbd);
  }
  else {
    printf("( A^T B^T)^T is not of full rank\n");
  }

END:
  nAG_FREE(a);
  nAG_FREE(alpha);
  nAG_FREE(b);
  nAG_FREE(beta);
  nAG_FREE(q);
  nAG_FREE(u);
  nAG_FREE(v);
  nAG_FREE(iwork);

  return exit_status;
}

#undef B
#undef A


関連情報
© 譌・譛ャ繝九Η繝シ繝。繝ェ繧ォ繝ォ繧「繝ォ繧エ繝ェ繧コ繝�繧コ繧ー繝ォ繝シ繝玲�ェ蠑丈シ夂、セ 2024
Privacy Policy  /  Trademarks