solve.cpp
Go to the documentation of this file.
1 #ifndef _USE_MATH_DEFINES
2 #define _USE_MATH_DEFINES
3 #endif
4 #include <cmath>
5 #include <GfxTL/MathHelper.h>
6
7 int dmat_solve(int n, int rhs_num, double a[])
8 //******************************************************************************
9 //
10 // Purpose:
11 //
12 // DMAT_SOLVE uses Gauss-Jordan elimination to solve an N by N linear system.
13 //
14 // Discussion:
15 //
16 // The doubly dimensioned array A is treated as a one dimensional vector,
17 // stored by COLUMNS. Entry A(I,J) is stored as A[I+J*N]
18 //
19 // Modified:
20 //
21 // 29 August 2003
22 //
23 // Author:
24 //
25 // John Burkardt
26 //
27 // Parameters:
28 //
29 // Input, int N, the order of the matrix.
30 //
31 // Input, int RHS_NUM, the number of right hand sides. RHS_NUM
32 // must be at least 0.
33 //
34 // Input/output, double A[N*(N+RHS_NUM)], contains in rows and columns 1
35 // to N the coefficient matrix, and in columns N+1 through
36 // N+RHS_NUM, the right hand sides. On output, the coefficient matrix
37 // area has been destroyed, while the right hand sides have
38 // been overwritten with the corresponding solutions.
39 //
40 // Output, int DMAT_SOLVE, singularity flag.
41 // 0, the matrix was not singular, the solutions were computed;
42 // J, factorization failed on step J, and the solutions could not
43 // be computed.
44 //
45 {
46  double apivot;
47  double factor;
48  int i;
49  int ipivot;
50  int j;
51  int k;
52  double temp;
53
54  for (j = 0; j < n; j++)
55  {
56  //
57  // Choose a pivot row.
58  //
59  ipivot = j;
60  apivot = a[j + j * n];
61
62  for (i = j; i < n; i++)
63  {
64  if (fabs(apivot) < fabs(a[i + j * n]))
65  {
66  apivot = a[i + j * n];
67  ipivot = i;
68  }
69  }
70
71  if (apivot == 0.0)
72  {
73  return j;
74  }
75  //
76  // Interchange.
77  //
78  for (i = 0; i < n + rhs_num; i++)
79  {
80  temp = a[ipivot + i * n];
81  a[ipivot + i * n] = a[j + i * n];
82  a[j + i * n] = temp;
83  }
84  //
85  // A(J,J) becomes 1.
86  //
87  a[j + j * n] = 1.0;
88  for (k = j; k < n + rhs_num; k++)
89  {
90  a[j + k * n] = a[j + k * n] / apivot;
91  }
92  //
93  // A(I,J) becomes 0.
94  //
95  for (i = 0; i < n; i++)
96  {
97  if (i != j)
98  {
99  factor = a[i + j * n];
100  a[i + j * n] = 0.0;
101  for (k = j; k < n + rhs_num; k++)
102  {
103  a[i + k * n] = a[i + k * n] - factor * a[j + k * n];
104  }
105  }
106
107  }
108
109  }
110
111  return 0;
112 }
113 //******************************************************************************
armarx::ctrlutil::a
double a(double t, double a0, double j)
Definition: CtrlUtil.h:45
dmat_solve
int dmat_solve(int n, int rhs_num, double a[])
Definition: solve.cpp:7
MathHelper.h