Math589A_GaussElimCAndPython/gauss_solve.c at main · Arizona-Math/Math589A_GaussElimCAndPython · GitHub

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
/*----------------------------------------------------------------
* File:     gauss_solve.c
*----------------------------------------------------------------
*
* Author:   Marek Rychlik (rychlik@arizona.edu)
* Date:     Sun Sep 22 15:40:29 2024
* Copying:  (C) Marek Rychlik, 2020. All rights reserved.
*
*----------------------------------------------------------------*/

#include <math.h>
#include <assert.h>
#include <stdio.h>
#include <string.h>

#include "gauss_solve.h"
#include "helpers.h"


void gauss_solve_in_place(const int n, double A[n][n], double b[n])
{
  for(int k = 0; k < n; ++k) {
    for(int i = k+1; i < n; ++i) {
      /* Store the multiplier into A[i][k] as it would become 0 and be
	 useless */
      A[i][k] /= A[k][k];
      for( int j = k+1; j < n; ++j) {
	A[i][j] -= A[i][k] * A[k][j];
      }
      b[i] -= A[i][k] * b[k];
    }
  } /* End of Gaussian elimination, start back-substitution. */
  for(int i = n-1; i >= 0; --i) {
    for(int j = i+1; j<n; ++j) {
      b[i] -= A[i][j] * b[j];
    }
    b[i] /= A[i][i];
  } /* End of back-substitution. */
}

void lu_in_place(const int n, double A[n][n])
{
  for(int k = 0; k < n; ++k) {
    for(int i = k; i < n; ++i) {
      for(int j=0; j<k; ++j) {
	/* U[k][i] -= L[k][j] * U[j][i] */
	A[k][i] -=  A[k][j] * A[j][i];
      }
    }
    for(int i = k+1; i<n; ++i) {
      for(int j=0; j<k; ++j) {
	/* L[i][k] -= A[i][k] * U[j][k] */
	A[i][k] -= A[i][j]*A[j][k];
      }
      /* L[i][k] /= U[k][k] */
      A[i][k] /= A[k][k];
    }
  }
}

/* Simply finds A = L.U from the packed form A = L+U */
void lu_in_place_reconstruct(int n, double A[n][n])
{
  for(int k = n-1; k >= 0; --k) {
    for(int i = k+1; i<n; ++i) {
      A[i][k] *= A[k][k];
      for(int j=0; j<k; ++j) {
	A[i][k] += A[i][j]*A[j][k];
      }
    }
    for(int i = k; i < n; ++i) {
      for(int j=0; j<k; ++j) {
	A[k][i] +=  A[k][j] * A[j][i];
      }
    }
  }
}

void plu_in_place(const int n, double A[n][n], int P[n])
{
  /* (Re)Initialize P */
  for(int k = 0; k < n; ++k) P[k] = k;

  for(int k = 0; k < n; ++k) {
    /* Find the pivot element */
    double pivot = A[k][k]; int r = k;

    for(int i = k + 1; i < n; ++i) {
      if( fabs(pivot) < fabs(A[i][k]) ) {
	pivot = A[i][k];
	r = i;
      }
    }

    /* row exchange */
    if(r != k)  {
      SWAP(P[k], P[r], int);

      /* Swap rows k and r in L, only cols 1 to k-1 */
      for(int j = 0; j < k; ++j) {
	// SWAP(L[k][j], L[r][j], double);
	SWAP(A[k][j], A[r][j], double);
      }
      /* In this version, we swap the rest of the row, too */
      for(int j = k; j < n; ++j) {
	// SWAP(U[k][j], U[r][j], double);
	SWAP(A[k][j], A[r][j], double);
      }
    }

    for(int i = k+1; i < n; ++i) {
      // L[i][k] = A[i][k] / A[k][k];
      A[i][k] /= A[k][k];
      for(int j = k+1; j < n; ++j) {
	// A[i][j] -=  L[i][k] * A[k][j];
	A[i][j] -=  A[i][k] * A[k][j];
      }
    }
  }
}