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
|
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
int rowcmp(float *row0, float *row1, int n)
{
for (int i = 0; i < n; ++i) {
if (row0[i])
return 0;
if (row1[i])
return 1;
}
return 0;
}
/* Upper-triangularise matrix mat of size mxn */
void ut(float *mat, int m, int n)
{
long i, j, k;
float *current, *other, *tmp = malloc(n * sizeof(float));
for (i = 0; i + 1 < m; ++i) {
current = mat + i * n;
for (j = i + 1; j < m; ++j) {
other = mat + j * n;
if (rowcmp(current, other, n)) {
memcpy(tmp, current, n);
memcpy(current, other, n);
memcpy(other, tmp, n);
}
}
for (j = i + 1; j < m; ++j) {
for (k = 0; k < n && !mat[i*n + k]; ++k);
if (k == n)
break;
float c = mat[j*n + k] / mat[i*n + k];
mat[j*n + k] = 0.0; /* floating point imprecision */
while (++k < n)
mat[j*n + k] -= c * mat[i*n + k];
}
}
free(tmp);
}
float *transpose(float *mat, int m, int n)
{
float *result = malloc(m * n * sizeof(float));
for (long i = 0; i < m; ++i)
for (long j = 0; j < n; ++j)
result[j*m + i] = mat[i*n + j];
return result;
}
void printmatf(char *format, float *mat, int m, int n)
{
for (long i = 0; i < m; ++i) {
for (long j = 0; j < n; ++j)
printf(format, mat[i*n + j]);
putchar(10);
}
}
int main()
{
int n;
puts("Please input the size of the square matrix:");
fputs("n = ", stdout);
scanf("%d", &n);
long i, j, k, len = n * n * sizeof(float);
float sum = 0.0, *a = malloc(len);
puts("Now please input the matrix:");
for (i = 0; i < n * n; ++i) {
scanf("%f", a + i);
sum += a[i];
}
printf("The sum of all the matrix elements: %.2f\n", sum);
float *b = transpose(a, n, n);
puts("The transpose of the matrix:");
printmatf("%8.3f", b, n, n);
free(a);
b = realloc(b, len * 2);
for (i = 0; i < n; ++i)
for (j = 0; j < n; ++j)
b[n*n + i*n + j] = i == j;
a = transpose(b, n * 2, n);
ut(a, n, n * 2);
free(b);
float det = 1.0;
for (i = 0; i < n; ++i)
det *= a[i*n*2 + i];
printf("The determinant of the matrix: %.2f\n", det);
if (!det) {
puts("The given matrix is not invertible.");
free(a);
return 0;
}
for (i = n - 1; i >= 0; --i) {
float d = a[i*n*2 + i];
for (j = n; j < n * 2; ++j)
a[i*n*2 + j] /= d;
for (k = i - 1; k >= 0; --k) {
float c = a[k*n*2 + i];
for (j = n; j < n * 2; ++j)
a[k*n*2 + j] -= a[i*n*2 + j] * c;
}
}
b = transpose(a, n, n * 2);
a = transpose(b + n * n, n, n);
puts("The inverse of the matrix:");
printmatf("%8.3f", a, n, n);
free(a);
free(b);
return 0;
}
|
