Keyword: 実行列, 誤差限界, 改良
概要
本サンプルは実連立一次方程式の誤差限界をもつ解の改良を行うC言語によるサンプルプログラムです。 本サンプルは行列Aと行列Bが以下で示される場合の連立一次方程式AX=Bについて反復改良を用いて解を改良し、算出された解と前進／後退誤差限界を出力します。
※本サンプルはnAG Cライブラリに含まれる関数 nag_dgerfs() のExampleコードです。本サンプル及び関数の詳細情報は nag_dgerfs のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_dgerfs のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10
| このデータをダウンロード |
nag_dgerfs (f07ahc) Example Program Data 4 2 :Values of N and NRHS 1.80 2.88 2.05 -0.89 5.25 -2.95 -0.95 -3.80 1.58 -2.69 -2.90 -1.04 -1.11 -0.66 -0.59 0.80 :End of matrix A 9.52 18.47 24.35 2.25 0.77 -13.28 -6.22 -6.21 :End of matrix B
- １行目はタイトル行で読み飛ばされます。
- ２行目に行列Aの次数(n)と右辺の数(nrhs)を指定しています。
- ３～６行目に行列Aの要素を指定しています。
- ７～１０行目に行列Bの要素を指定しています。
出力結果
（本関数の詳細はnag_dgerfs のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10 11 12 13 14
| この出力例をダウンロード |
nag_dgerfs (f07ahc) Example Program Results
Solution(s)
1 2
1 1.0000 3.0000
2 -1.0000 2.0000
3 3.0000 4.0000
4 -5.0000 1.0000
Backward errors (machine-dependent)
5.6e-17 6.2e-17
Estimated forward error bounds (machine-dependent)
2.4e-14 3.3e-14
- ７～１０行目にxの解が出力されています。
- １３行目に後退誤差限界が出力されています。
- １５行目に推定された前進誤差限界が出力されています。
ソースコード
（本関数の詳細はnag_dgerfs のマニュアルページを参照）
※本サンプルソースコードは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
| このソースコードをダウンロード |
/* nag_dgerfs (f07ahc) 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 <nagx04.h>
int main(void)
{
/* Scalars */
Integer berr_len, i, ipiv_len, ferr_len, j, n, nrhs;
Integer pda, pdaf, pdb, pdx;
Integer exit_status = 0;
NagError fail;
Nag_OrderType order;
/* Arrays */
double *a = 0, *af = 0, *b = 0, *berr = 0, *ferr = 0, *x = 0;
Integer *ipiv = 0;
#ifdef nAG_COLUMN_MAJOR
#define A(I, J) a[(J-1)*pda + I - 1]
#define AF(I, J) af[(J-1)*pdaf + I - 1]
#define B(I, J) b[(J-1)*pdb + I - 1]
#define X(I, J) x[(J-1)*pdx + I - 1]
order = Nag_ColMajor;
#else
#define A(I, J) a[(I-1)*pda + J - 1]
#define AF(I, J) af[(I-1)*pdaf + J - 1]
#define B(I, J) b[(I-1)*pdb + J - 1]
#define X(I, J) x[(I-1)*pdx + J - 1]
order = Nag_RowMajor;
#endif
INIT_FAIL(fail);
printf("nag_dgerfs (f07ahc) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%ld%ld%*[^\n] ", &n, &nrhs);
#ifdef nAG_COLUMN_MAJOR
pda = n;
pdaf = n;
pdb = n;
pdx = n;
#else
pda = n;
pdaf = n;
pdb = nrhs;
pdx = nrhs;
#endif
berr_len = nrhs;
ferr_len = nrhs;
ipiv_len = n;
/* Allocate memory */
if (!(a = nAG_ALLOC(n * n, double)) ||
!(af = nAG_ALLOC(n * n, double)) ||
!(b = nAG_ALLOC(n * nrhs, double)) ||
!(berr = nAG_ALLOC(berr_len, double)) ||
!(ferr = nAG_ALLOC(ferr_len, double)) ||
!(x = nAG_ALLOC(n * nrhs, double)) ||
!(ipiv = nAG_ALLOC(ipiv_len, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Read A and B from data file, and copy A to AF and B to X */
for (i = 1; i <= n; ++i) {
for (j = 1; j <= n; ++j)
scanf("%lf", &A(i, j));
}
scanf("%*[^\n] ");
for (i = 1; i <= n; ++i) {
for (j = 1; j <= nrhs; ++j)
scanf("%lf", &B(i, j));
}
scanf("%*[^\n] ");
for (i = 1; i <= n; ++i) {
for (j = 1; j <= n; ++j)
AF(i, j) = A(i, j);
}
for (i = 1; i <= n; ++i) {
for (j = 1; j <= nrhs; ++j)
X(i, j) = B(i, j);
}
/* Factorize A in the array AF */
/* nag_dgetrf (f07adc).
* LU factorization of real m by n matrix
*/
nag_dgetrf(order, n, n, af, pdaf, ipiv, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_dgetrf (f07adc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf("\n");
/* Compute solution in the array X */
/* nag_dgetrs (f07aec).
* Solution of real system of linear equations, multiple
* right-hand sides, matrix already factorized by nag_dgetrf
* (f07adc)
*/
nag_dgetrs(order, Nag_NoTrans, n, nrhs, af, pdaf, ipiv, x, pdx, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_dgetrs (f07aec).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Improve solution, and compute backward errors and */
/* estimated bounds on the forward errors */
/* nag_dgerfs (f07ahc).
* Refined solution with error bounds of real system of
* linear equations, multiple right-hand sides
*/
nag_dgerfs(order, Nag_NoTrans, n, nrhs, a, pda, af, pdaf, ipiv, b, pdb, x,
pdx, ferr, berr, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_dgerfs (f07ahc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Print solution */
/* 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,
x, pdx, "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;
}
printf("\nBackward errors (machine-dependent)\n");
for (j = 1; j <= nrhs; ++j)
printf("%11.1e%s", berr[j - 1], j % 7 == 0 ? "\n" : " ");
printf("\nEstimated forward error bounds (machine-dependent)\n");
for (j = 1; j <= nrhs; ++j)
printf("%11.1e%s", ferr[j - 1], j % 7 == 0 ? "\n" : " ");
printf("\n");
END:
nAG_FREE(a);
nAG_FREE(af);
nAG_FREE(b);
nAG_FREE(berr);
nAG_FREE(ferr);
nAG_FREE(x);
nAG_FREE(ipiv);
return exit_status;
}