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