Matrix-Algebra

1
#include <iostream>
2
#include <iomanip>
3
#include <cmath>
4

5
using namespace std;
6

7
class MatrixException : public exception {
8
private:
9
    char *message;
10

11
public:
12
    MatrixException(char *msg) : message(msg) {}
13

14
    char *what() {
15
        return message;
16
    }
17
};
18

19

20
class Matrix {
21
protected:
22
    unsigned int n;
23
    unsigned int m;
24
    double **data;
25

26
public:
27
    Matrix(int cols, int rows) {
28
        n = cols;
29
        m = rows;
30

31
        data = new double *[cols];
32
        for (int i = 0; i < cols; i++) {
33
            data[i] = new double[rows];
34
            for (int j = 0; j < rows; j++) {
35
                data[i][j] = 0;
36
            }
37
        }
38
    }
39

40

41
    int getN() {
42
        return n;
43
    }
44

45
    int getM() {
46
        return m;
47
    }
48

49

50
    void setValue(int i, int j, double value) {
51
        data[i][j] = value;
52
    }
53

54
    double getValue(int i, int j) {
55
        return data[i][j];
56
    }
57

58
//    double *getRow(int i) {
59
//        double ans[n];
60
//        for (int k = 0; k < n; k++) {
61
//            ans[k] = data[i][k];
62
//        }
63
//
64
//        return ans;
65
//    }
66

67
    Matrix getCopy() {
68
        Matrix ans = Matrix(getN(), getM());
69
        for (int i = 0; i < getN(); i++) {
70
            for (int j = 0; j < getM(); j++) {
71
                ans.setValue(i, j, getValue(i, j));
72
            }
73
        }
74
        return ans;
75
    }
76

77
    Matrix operator+(Matrix const secondMatrix) const {
78
        if (n != secondMatrix.n || m != secondMatrix.m) {
79
            throw MatrixException((char *) "Error: the dimensional problem occurred");
80
        }
81
        Matrix answer = Matrix(n, m);
82
        for (int i = 0; i < n; i++) {
83
            for (int j = 0; j < m; j++) {
84
                answer.data[i][j] = data[i][j] + secondMatrix.data[i][j];
85
            }
86
        }
87
        return answer;
88
    }
89

90
    Matrix operator-(Matrix const secondMatrix) const {
91
        if (n != secondMatrix.n || m != secondMatrix.m) {
92
            throw MatrixException((char *) "Error: the dimensional problem occurred");
93
        }
94
        Matrix answer = Matrix(n, m);
95
        for (int i = 0; i < n; i++) {
96
            for (int j = 0; j < m; j++) {
97
                answer.data[i][j] = data[i][j] - secondMatrix.data[i][j];
98
            }
99
        }
100
        return answer;
101
    }
102

103
    Matrix operator*(Matrix const secondMatrix) const {
104
        if (m != secondMatrix.n) {
105
            throw MatrixException((char *) "Error: the dimensional problem occurred");
106
        }
107
        Matrix answer = Matrix(n, secondMatrix.m);
108

109
        for (int i = 0; i < n; i++) {
110
            for (int j = 0; j < secondMatrix.m; j++) {
111

112
                for (int k = 0; k < secondMatrix.n; k++) {
113
                    answer.setValue(i, j, answer.data[i][j] + data[i][k] * secondMatrix.data[k][j]);
114
                }
115
            }
116

117
        }
118

119
        return answer;
120
    }
121

122
    Matrix getTranspose() {
123

124
        Matrix answer = Matrix(m, n);
125

126
        for (int i = 0; i < n; i++) {
127
            for (int j = 0; j < m; j++) {
128
                answer.setValue(j, i, data[i][j]);
129
            }
130
        }
131

132
        return answer;
133

134
    }
135

136
    Matrix &operator=(Matrix const &otherMatrix) {
137
        data = otherMatrix.data;
138
        n = otherMatrix.n;
139
        m = otherMatrix.m;
140
        return *this;
141
    }
142

143
    friend istream &operator>>(istream &is, Matrix &matrix) {
144
        for (int i = 0; i < matrix.n; i++) {
145
            for (int j = 0; j < matrix.m; j++) {
146
                double t;
147
                is >> t;
148
                matrix.setValue(i, j, t);
149
            }
150
        }
151
        return is;
152
    }
153

154
    friend ostream &operator<<(ostream &os, Matrix &matrix) {
155
        for (int i = 0; i < matrix.n; i++) {
156
            for (int j = 0; j < matrix.m; j++) {
157
                double num = matrix.getValue(i, j);
158
                if (abs(num) < 0.005) {
159
                    os << setprecision(4) << fixed << abs(matrix.getValue(i, j));
160
                } else {
161
                    os << setprecision(4) << fixed << matrix.getValue(i, j);
162
                }
163
                if (j != matrix.m - 1) {
164
                    cout << " ";
165
                }
166

167

168
            }
169
            os << endl;
170
        }
171
        return os;
172
    }
173

174
    double getDeterminant();
175

176
    Matrix getInverse();
177

178
    Matrix solveGaussianEquation(Matrix b);
179

180
    Matrix solveJacobianEquation(Matrix b, long double epsilon);
181

182
    Matrix solveSeidelEquation(Matrix b, long double epsilon);
183

184

185
    static void printAugmented(Matrix left, Matrix right) {
186

187
        for (int i = 0; i < left.getN(); i++) {
188
            for (int j = 0; j < left.getM(); j++) {
189
                cout << setprecision(2) << fixed << left.getValue(i, j) << " ";
190
            }
191
            for (int j = 0; j < right.getM(); j++) {
192
                cout << setprecision(2) << fixed << right.getValue(i, j) << " ";
193
            }
194
            cout << endl;
195
        }
196
    }
197

198
};
199

200
class ColumnVector : public Matrix {
201
public:
202
    ColumnVector(int n) : Matrix(n, 1) {
203

204
    }
205

206
};
207

208
class SquareMatrix : public Matrix {
209
public:
210
    SquareMatrix(int n) : Matrix(n, n) {
211

212
    }
213
};
214

215
class IdentityMatrix : public SquareMatrix {
216
public:
217
    IdentityMatrix(int n) : SquareMatrix(n) {
218
        for (int i = 0; i < n; i++) {
219
            for (int j = 0; j < n; j++) {
220
                if (i == j) {
221
                    setValue(i, j, 1);
222
                }
223
            }
224
        }
225
    }
226

227
    IdentityMatrix(Matrix matrix) : SquareMatrix(max(matrix.getN(), matrix.getM())) {
228
        int maxSize = max(matrix.getN(), matrix.getM());
229
        int minSize = min(matrix.getN(), matrix.getM());
230
        for (int i = 0; i < maxSize; i++) {
231
            for (int j = 0; j < maxSize; j++) {
232
                if (i == j && i < minSize) {
233
                    setValue(i, j, 1);
234
                }
235
            }
236
        }
237
    }
238
};
239

240
class EliminationMatrix : public IdentityMatrix {
241
public:
242
    EliminationMatrix(Matrix matrix, int i, int j) : IdentityMatrix(matrix.getN()) {
243
        i--;
244
        j--;
245
        double e = -matrix.getValue(i, j) / matrix.getValue(j, j);
246
        setValue(i, j, e);
247
    }
248
};
249

250
class PermutationMatrix : public IdentityMatrix {
251
public:
252
    PermutationMatrix(Matrix matrix, int i, int j) : IdentityMatrix(matrix.getN()) {
253
        i--;
254
        j--;
255
        setValue(i, i, 0);
256
        setValue(j, j, 0);
257
        setValue(i, j, 1);
258
        setValue(j, i, 1);
259
    }
260
};
261

262
double Matrix::getDeterminant() {
263
    Matrix result = getCopy();
264
    int steps = 1;
265

266
    for (int j = 0; j < getM(); j++) {
267
        if (j != result.getM() - 1) {
268
            double pivot = result.getValue(j, j);
269
            int pivot_i = j;
270
            for (int i = j; i < getN(); i++) {
271
                if (abs(result.getValue(i, j)) > abs(pivot)) {
272
                    pivot = result.getValue(i, j);
273
                    pivot_i = i;
274
                }
275
            }
276
            if (pivot_i != j) {
277
                Matrix p = PermutationMatrix(result, pivot_i + 1, j + 1);
278
                result = p * result;
279
                cout << "step #" << steps << ": permutation" << endl;
280
                cout << result;
281
                steps++;
282
            }
283
        }
284
        for (int i = j + 1; i < getN(); i++) {
285
            if (result.getValue(i, j) != 0) {
286
                Matrix e = EliminationMatrix(result, i + 1, j + 1);
287
                result = e * result;
288
                cout << "step #" << steps << ": elimination" << endl;
289
                cout << result;
290
                steps++;
291
            }
292
        }
293
    }
294

295
    double det = 1;
296
    for (int i = 0; i < n; i++) {
297
        det *= result.getValue(i, i);
298
    }
299
    cout << "result:" << endl << setprecision(2) << fixed << det << endl;
300
    return det;
301
}
302

303
Matrix Matrix::getInverse() {
304
    if (getN() != getM()) {
305
        throw MatrixException((char *) "Error: the dimensional problem occurred");
306
    }
307
    Matrix left = getCopy();
308
    Matrix right = IdentityMatrix(getN());
309

310
    for (int j = 0; j < getM(); j++) {
311
        if (j != left.getM() - 1) {
312
            double pivot = left.getValue(j, j);
313
            int pivot_i = j;
314
            for (int i = j; i < getN(); i++) {
315
                if (abs(left.getValue(i, j)) > abs(pivot)) {
316
                    pivot = left.getValue(i, j);
317
                    pivot_i = i;
318
                }
319
            }
320
            if (pivot_i != j) {
321
                Matrix p = PermutationMatrix(left, pivot_i + 1, j + 1);
322
                left = p * left;
323
                right = p * right;
324
            }
325
        }
326
        for (int i = j + 1; i < getN(); i++) {
327
            if (left.getValue(i, j) != 0) {
328
                Matrix e = EliminationMatrix(left, i + 1, j + 1);
329
                left = e * left;
330
                right = e * right;
331
            }
332
        }
333
    }
334

335
    for (int j = getM() - 1; j >= 0; j--) {
336

337
        for (int i = j - 1; i >= 0; i--) {
338
            if (left.getValue(i, j) != 0) {
339
                Matrix e = EliminationMatrix(left, i + 1, j + 1);
340
                left = e * left;
341
                right = e * right;
342
            }
343
        }
344
    }
345

346
    for (int i = 0; i < getN(); i++) {
347
        double value = 1 / left.getValue(i, i);
348
        left.setValue(i, i, left.getValue(i, i) * value);
349

350
        for (int j = 0; j < getN(); j++) {
351
            right.setValue(i, j, right.getValue(i, j) * value);
352
        }
353
    }
354

355

356
    return right;
357
}
358

359
Matrix Matrix::solveGaussianEquation(Matrix b) {
360
    if (b.getM() != 1) {
361
        if (getN() != getM()) {
362
            throw MatrixException((char *) "Error: the dimensional problem occurred");
363
        }
364
    }
365
    Matrix left = getCopy();
366
    Matrix right = b.getCopy();
367

368

369
    for (int j = 0; j < getM(); j++) {
370
        if (j != left.getM() - 1) {
371
            double pivot = left.getValue(j, j);
372
            int pivot_i = j;
373
            for (int i = j; i < getN(); i++) {
374
                if (abs(left.getValue(i, j)) > abs(pivot)) {
375
                    pivot = left.getValue(i, j);
376
                    pivot_i = i;
377
                }
378
            }
379
            if (pivot_i != j) {
380
                Matrix p = PermutationMatrix(left, pivot_i + 1, j + 1);
381
                left = p * left;
382
                right = p * right;
383
            }
384
        }
385
        for (int i = j + 1; i < getN(); i++) {
386
            if (left.getValue(i, j) != 0) {
387
                Matrix e = EliminationMatrix(left, i + 1, j + 1);
388
                left = e * left;
389
                right = e * right;
390
            }
391
        }
392
    }
393
    for (int j = getM() - 1; j >= 0; j--) {
394
        for (int i = j - 1; i >= 0; i--) {
395
            if (left.getValue(i, j) != 0) {
396
                Matrix e = EliminationMatrix(left, i + 1, j + 1);
397
                left = e * left;
398
                right = e * right;
399
            }
400
        }
401
    }
402

403
    for (int i = 0; i < getN(); i++) {
404
        double value = 1 / left.getValue(i, i);
405

406
        left.setValue(i, i, left.getValue(i, i) * value);
407
        right.setValue(i, 0, right.getValue(i, 0) * value);
408
    }
409

410
    return right;
411
}
412

413
Matrix Matrix::solveJacobianEquation(Matrix b, long double epsilon) {
414
    if (getM() != b.getN()) {
415
        throw MatrixException("The method is not applicable!");
416
    }
417
    for (int i = 0; i < getN(); i++) {
418
        if (getValue(i, i) == 0) {
419
            throw MatrixException("The method is not applicable!");
420
        }
421
        long double sum = 0;
422
        for (int j = 0; j < getM(); j++) {
423
            if (i != j) {
424
                sum += getValue(i, j);
425
            }
426
        }
427
        if (getValue(i, i) < sum) {
428
            throw MatrixException("The method is not applicable!");
429
        }
430
    }
431
    Matrix aD = SquareMatrix(n);
432
    for (int i = 0; i < getN(); i++) {
433
        aD.setValue(i, i, getValue(i, i));
434
    }
435
    Matrix aDInverse = aD.getInverse();
436

437
    Matrix alpha = aDInverse * (aD - *this);
438
    cout << "alpha:" << endl << alpha;
439

440
    Matrix beta = aDInverse * b;
441
    cout << "beta:" << endl << beta;
442

443
    Matrix x = beta.getCopy();
444
    Matrix newx = IdentityMatrix(beta);
445
    cout << "x(0):" << endl << x;
446
    int steps = 1;
447

448
    while (true) {
449
        newx = alpha * x + beta;
450
        Matrix delta = newx - x;
451
        long double error = 0;
452
        for (int i = 0; i < delta.getN(); i++) {
453
            error += powl(delta.getValue(i, 0), 2);
454
        }
455
        error = ::sqrt(error);
456
        cout << "e: " << error << endl;
457

458

459
        cout << "x(" << steps++ << "):" << endl << newx;
460

461
        if (error <= epsilon) {
462
            return x;
463
        }
464

465
        x = newx;
466
    }
467
}
468

469
Matrix Matrix::solveSeidelEquation(Matrix b, long double epsilon) {
470
    if (getM() != b.getN()) {
471
        throw MatrixException("The method is not applicable!");
472
    }
473
    for (int i = 0; i < getN(); i++) {
474
        if (getValue(i, i) == 0) {
475
            throw MatrixException("The method is not applicable!");
476
        }
477
        long double sum = 0;
478
        for (int j = 0; j < getM(); j++) {
479
            if (i != j) {
480
                sum += fabs(getValue(i, j));
481
            }
482
        }
483
        if (fabs(getValue(i, i)) < sum) {
484
            throw MatrixException("The method is not applicable!");
485
        }
486
    }
487

488
    Matrix aD = SquareMatrix(n);
489
    for (int i = 0; i < getN(); i++) {
490
        aD.setValue(i, i, getValue(i, i));
491
    }
492

493
    Matrix aDInverse = aD.getInverse();
494

495
    Matrix alpha = aDInverse * (aD - *this);
496
    Matrix beta = aDInverse * b;
497

498
    Matrix aL = SquareMatrix(n);
499
    for (int j = 0; j < getM(); j++){
500
        for (int i = j + 1; i < getN(); i++){
501
            aL.setValue(i, j, alpha.getValue(i, j));
502
        }
503
    }
504

505
    Matrix aU = SquareMatrix(n);
506
    for (int j = getM() - 1; j >= 0; j--) {
507
        for (int i = j - 1; i >= 0; i--){
508
            aU.setValue(i, j, alpha.getValue(i, j));
509
        }
510
    }
511

512
    Matrix IminusB = IdentityMatrix(aL) - aL;
513
    Matrix IminusBInverse = IminusB.getInverse();
514

515
    cout << "beta:" << endl << beta;
516
    cout << "alpha:" << endl << alpha;
517
    cout << "B:" << endl << aL;
518
    cout << "C:" << endl << aU;
519
    cout << "I-B:" << endl << IminusB;
520
    cout << "(I-B)_-1:" << endl << IminusBInverse;
521
    cout << "x(0):" << endl << beta;
522

523
    Matrix x = beta.getCopy();
524
    int steps = 1;
525
    while (true){
526
        Matrix newx = IminusBInverse * (aU * x + beta);
527

528
        Matrix delta = newx - x;
529
        long double error = 0;
530
        for (int i = 0; i < delta.getN(); i++) {
531
            error += powl(delta.getValue(i, 0), 2);
532
        }
533
        error = ::sqrt(error);
534
        cout << "e: " << error << endl;
535

536
        cout << "x(" << steps++ << "):" << endl << newx;
537

538
        if (error <= epsilon) {
539
            return x;
540
        }
541

542
        x = newx;
543
    }
544
}
545

546
int main() {
547
    int n;
548
    cin >> n;
549

550
    Matrix a = SquareMatrix(n);
551
    cin >> a;
552

553
    cin >> n;
554
    Matrix b = ColumnVector(n);
555
    cin >> b;
556

557
    double epsilon;
558
    cin >> epsilon;
559

560
    try {
561
        a.solveSeidelEquation(b, epsilon);
562
    }
563
    catch (MatrixException e) {
564
        cout << e.what() << endl;
565
    }
566
    return 0;
567

568
}
569

570

571

572