Orbiter 2022
Combinatorial Objects
linear_algebra_RREF.cpp
Go to the documentation of this file.
1/*
2 * linear_algebra_RREF.cpp
3 *
4 * Created on: Feb 7, 2021
5 * Author: betten
6 */
7
8
9
10#include "foundations.h"
11
12using namespace std;
13
14
15namespace orbiter {
16namespace layer1_foundations {
17namespace linear_algebra {
18
19
21 int f_special, int f_complete, int *base_cols,
22 int f_P, int *P, int m, int n, int Pn, int verbose_level)
23// returns the rank which is the number of entries in base_cols
24// A is a m x n matrix,
25// P is a m x Pn matrix (if f_P is TRUE)
26{
27 int f_v = (verbose_level >= 1);
28 int f_vv = FALSE; //(verbose_level >= 2);
29 int f_vvv = FALSE; //(verbose_level >= 3);
30 int rank, i, j, k, jj;
31 int pivot, pivot_inv = 0, a, b, c, z, f;
33
34 if (f_v) {
35 cout << "linear_algebra::Gauss_int" << endl;
36 }
37 if (f_vv) {
38 cout << "Gauss algorithm for matrix:" << endl;
39 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
40 //print_tables();
41 }
42 i = 0;
43 for (j = 0; j < n; j++) {
44 if (f_vv) {
45 cout << "searching for pivot element in column j=" << j << endl;
46 }
47 // search for pivot element:
48 for (k = i; k < m; k++) {
49 if (A[k * n + j]) {
50 if (f_vv) {
51 cout << "i=" << i << " pivot found in "
52 << k << "," << j << endl;
53 }
54 // pivot element found:
55 if (k != i) {
56 for (jj = j; jj < n; jj++) {
57 Algo.int_swap(A[i * n + jj], A[k * n + jj]);
58 }
59 if (f_P) {
60 for (jj = 0; jj < Pn; jj++) {
61 Algo.int_swap(P[i * Pn + jj], P[k * Pn + jj]);
62 }
63 }
64 }
65 break;
66 } // if != 0
67 } // next k
68
69 if (k == m) { // no pivot found
70 if (f_vv) {
71 cout << "no pivot found" << endl;
72 }
73 continue; // increase j, leave i constant
74 }
75
76 if (f_vv) {
77 cout << "row " << i << " pivot in row "
78 << k << " colum " << j << endl;
79 }
80
81 base_cols[i] = j;
82 //if (FALSE) {
83 // cout << ".";
84 // }
85
86 pivot = A[i * n + j];
87 if (f_vv) {
88 cout << "pivot=" << pivot << endl;
89 }
90 //pivot_inv = inv_table[pivot];
91 pivot_inv = F->inverse(pivot);
92 if (f_vv) {
93 cout << "pivot=" << pivot << " pivot_inv="
94 << pivot_inv << endl;
95 }
96 if (!f_special) {
97 // make pivot to 1:
98 for (jj = j; jj < n; jj++) {
99 A[i * n + jj] = F->mult(A[i * n + jj], pivot_inv);
100 }
101 if (f_P) {
102 for (jj = 0; jj < Pn; jj++) {
103 P[i * Pn + jj] = F->mult(P[i * Pn + jj], pivot_inv);
104 }
105 }
106 if (f_vv) {
107 cout << "pivot=" << pivot << " pivot_inv=" << pivot_inv
108 << " made to one: " << A[i * n + j] << endl;
109 }
110 if (f_vvv) {
111 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
112 }
113 }
114
115 // do the gaussian elimination:
116
117 if (f_vv) {
118 cout << "doing elimination in column " << j << " from row "
119 << i + 1 << " to row " << m - 1 << ":" << endl;
120 }
121 for (k = i + 1; k < m; k++) {
122 if (f_vv) {
123 cout << "k=" << k << endl;
124 }
125 z = A[k * n + j];
126 if (z == 0) {
127 continue;
128 }
129 if (f_special) {
130 f = F->mult(z, pivot_inv);
131 }
132 else {
133 f = z;
134 }
135 f = F->negate(f);
136 A[k * n + j] = 0;
137 if (f_vv) {
138 cout << "eliminating row " << k << endl;
139 }
140 for (jj = j + 1; jj < n; jj++) {
141 a = A[i * n + jj];
142 b = A[k * n + jj];
143 // c := b + f * a
144 // = b - z * a if !f_special
145 // b - z * pivot_inv * a if f_special
146 c = F->mult(f, a);
147 c = F->add(c, b);
148 A[k * n + jj] = c;
149 if (f_vv) {
150 cout << A[k * n + jj] << " ";
151 }
152 }
153 if (f_P) {
154 for (jj = 0; jj < Pn; jj++) {
155 a = P[i * Pn + jj];
156 b = P[k * Pn + jj];
157 // c := b - z * a
158 c = F->mult(f, a);
159 c = F->add(c, b);
160 P[k * Pn + jj] = c;
161 }
162 }
163 if (f_vv) {
164 cout << endl;
165 }
166 if (f_vvv) {
167 cout << "A=" << endl;
168 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
169 }
170 }
171 i++;
172 if (f_vv) {
173 cout << "A=" << endl;
174 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
175 //print_integer_matrix(cout, A, m, n);
176 if (f_P) {
177 cout << "P=" << endl;
178 Int_vec_print_integer_matrix(cout, P, m, Pn);
179 }
180 }
181 } // next j
182 rank = i;
183 if (f_complete) {
184 //if (FALSE) {
185 // cout << ";"; cout.flush();
186 // }
187 for (i = rank - 1; i >= 0; i--) {
188 if (f_v) {
189 cout << ".";
190 }
191 j = base_cols[i];
192 if (!f_special) {
193 a = A[i * n + j];
194 }
195 else {
196 pivot = A[i * n + j];
197 pivot_inv = F->inverse(pivot);
198 }
199 // do the gaussian elimination in the upper part:
200 for (k = i - 1; k >= 0; k--) {
201 z = A[k * n + j];
202 if (z == 0) {
203 continue;
204 }
205 A[k * n + j] = 0;
206 for (jj = j + 1; jj < n; jj++) {
207 a = A[i * n + jj];
208 b = A[k * n + jj];
209 if (f_special) {
210 a = F->mult(a, pivot_inv);
211 }
212 c = F->mult(z, a);
213 c = F->negate(c);
214 c = F->add(c, b);
215 A[k * n + jj] = c;
216 }
217 if (f_P) {
218 for (jj = 0; jj < Pn; jj++) {
219 a = P[i * Pn + jj];
220 b = P[k * Pn + jj];
221 if (f_special) {
222 a = F->mult(a, pivot_inv);
223 }
224 c = F->mult(z, a);
225 c = F->negate(c);
226 c = F->add(c, b);
227 P[k * Pn + jj] = c;
228 }
229 }
230 } // next k
231 } // next i
232 }
233 if (f_vv) {
234 cout << endl;
235 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
236 //print_integer_matrix(cout, A, rank, n);
237 cout << "the rank is " << rank << endl;
238 }
239 if (f_v) {
240 cout << "linear_algebra::Gauss_int done" << endl;
241 }
242 return rank;
243}
244
245int linear_algebra::Gauss_int_with_pivot_strategy(int *A,
246 int f_special, int f_complete, int *pivot_perm,
247 int m, int n,
248 int (*find_pivot_function)(int *A, int m, int n, int r,
249 int *pivot_perm, void *data),
250 void *find_pivot_data,
251 int verbose_level)
252// returns the rank which is the number of entries in pivots
253// A is a m x n matrix
254{
255 int f_v = (verbose_level >= 1);
256 int f_vv = FALSE; //(verbose_level >= 2);
257 int f_vvv = FALSE; //(verbose_level >= 3);
258 int rank, i, j, k, jj;
259 int pivot, pivot_inv = 0, a, b, c, z, f, pi;
261
262 if (f_v) {
263 cout << "linear_algebra::Gauss_int_with_pivot_strategy" << endl;
264 }
265 if (f_vv) {
266 cout << "linear_algebra::Gauss_int_with_pivot_strategy "
267 "Gauss algorithm for matrix:" << endl;
268 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
269 //print_tables();
270 }
271 for (i = 0; i < m; i++) {
272 if (f_vv) {
273 cout << "i=" << i << endl;
274 }
275
276 j = (*find_pivot_function)(A, m, n, i, pivot_perm, find_pivot_data);
277
278 if (j == -1) {
279 break;
280 }
281
282 pi = pivot_perm[i];
283 pivot_perm[i] = j;
284 pivot_perm[j] = pi;
285
286 // search for pivot element in column j from row i down:
287 for (k = i; k < m; k++) {
288 if (A[k * n + j]) {
289 break;
290 } // if != 0
291 } // next k
292
293 if (k == m) { // no pivot found
294 if (f_vv) {
295 cout << "linear_algebra::Gauss_int_with_pivot_strategy "
296 "no pivot found in column " << j << endl;
297 }
298 exit(1);
299 }
300
301 if (f_vv) {
302 cout << "row " << i << " pivot in row " << k
303 << " colum " << j << endl;
304 }
305
306
307
308
309
310 // pivot element found in row k, check if we need to swap rows:
311 if (k != i) {
312 for (jj = 0; jj < n; jj++) {
313 Algo.int_swap(A[i * n + jj], A[k * n + jj]);
314 }
315 }
316
317
318 // now, pivot is in row i, column j :
319
320 pivot = A[i * n + j];
321 if (f_vv) {
322 cout << "pivot=" << pivot << endl;
323 }
324 pivot_inv = F->inverse(pivot);
325 if (f_vv) {
326 cout << "pivot=" << pivot << " pivot_inv=" << pivot_inv << endl;
327 }
328 if (!f_special) {
329 // make pivot to 1:
330 for (jj = 0; jj < n; jj++) {
331 A[i * n + jj] = F->mult(A[i * n + jj], pivot_inv);
332 }
333 if (f_vv) {
334 cout << "pivot=" << pivot << " pivot_inv=" << pivot_inv
335 << " made to one: " << A[i * n + j] << endl;
336 }
337 if (f_vvv) {
338 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
339 }
340 }
341
342 // do the gaussian elimination:
343
344 if (f_vv) {
345 cout << "doing elimination in column " << j << " from row "
346 << i + 1 << " down to row " << m - 1 << ":" << endl;
347 }
348 for (k = i + 1; k < m; k++) {
349 if (f_vv) {
350 cout << "k=" << k << endl;
351 }
352 z = A[k * n + j];
353 if (z == 0) {
354 continue;
355 }
356 if (f_special) {
357 f = F->mult(z, pivot_inv);
358 }
359 else {
360 f = z;
361 }
362 f = F->negate(f);
363 //A[k * n + j] = 0;
364 if (f_vv) {
365 cout << "eliminating row " << k << endl;
366 }
367 for (jj = 0; jj < n; jj++) {
368 a = A[i * n + jj];
369 b = A[k * n + jj];
370 // c := b + f * a
371 // = b - z * a if !f_special
372 // b - z * pivot_inv * a if f_special
373 c = F->mult(f, a);
374 c = F->add(c, b);
375 A[k * n + jj] = c;
376 if (f_vv) {
377 cout << A[k * n + jj] << " ";
378 }
379 }
380 if (f_vv) {
381 cout << endl;
382 }
383 if (f_vvv) {
384 cout << "A=" << endl;
385 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
386 }
387 }
388 i++;
389 if (f_vv) {
390 cout << "A=" << endl;
391 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
392 //print_integer_matrix(cout, A, m, n);
393 }
394 } // next j
395 rank = i;
396 if (f_complete) {
397 for (i = rank - 1; i >= 0; i--) {
398 j = pivot_perm[i];
399 if (!f_special) {
400 a = A[i * n + j];
401 }
402 else {
403 pivot = A[i * n + j];
404 pivot_inv = F->inverse(pivot);
405 }
406 // do the gaussian elimination in the upper part:
407 for (k = i - 1; k >= 0; k--) {
408 z = A[k * n + j];
409 if (z == 0) {
410 continue;
411 }
412 //A[k * n + j] = 0;
413 for (jj = 0; jj < n; jj++) {
414 a = A[i * n + jj];
415 b = A[k * n + jj];
416 if (f_special) {
417 a = F->mult(a, pivot_inv);
418 }
419 c = F->mult(z, a);
420 c = F->negate(c);
421 c = F->add(c, b);
422 A[k * n + jj] = c;
423 }
424 } // next k
425 } // next i
426 }
427 if (f_vv) {
428 cout << endl;
429 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
430 //print_integer_matrix(cout, A, rank, n);
431 cout << "the rank is " << rank << endl;
432 }
433 if (f_v) {
434 cout << "linear_algebra::Gauss_int_with_pivot_strategy done" << endl;
435 }
436 return rank;
437}
438
439int linear_algebra::Gauss_int_with_given_pivots(int *A,
440 int f_special, int f_complete, int *pivots, int nb_pivots,
441 int m, int n,
442 int verbose_level)
443// A is a m x n matrix
444{
445 int f_v = (verbose_level >= 1);
446 int f_vv = FALSE; //(verbose_level >= 2);
447 int f_vvv = FALSE; //(verbose_level >= 3);
448 int i, j, k, jj;
449 int pivot, pivot_inv = 0, a, b, c, z, f;
451
452 if (f_v) {
453 cout << "linear_algebra::Gauss_int_with_given_pivots" << endl;
454 }
455 if (f_vv) {
456 cout << "linear_algebra::Gauss_int_with_given_pivots "
457 "Gauss algorithm for matrix:" << endl;
458 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
459 cout << "pivots: ";
460 Int_vec_print(cout, pivots, nb_pivots);
461 cout << endl;
462 //print_tables();
463 }
464 for (i = 0; i < nb_pivots; i++) {
465 if (f_vv) {
466 cout << "i=" << i << endl;
467 }
468
469 j = pivots[i];
470
471 // search for pivot element in column j from row i down:
472 for (k = i; k < m; k++) {
473 if (A[k * n + j]) {
474 break;
475 } // if != 0
476 } // next k
477
478 if (k == m) { // no pivot found
479 if (f_v) {
480 cout << "linear_algebra::Gauss_int_with_given_pivots "
481 "no pivot found in column " << j << endl;
482 }
483 return FALSE;
484 }
485
486 if (f_vv) {
487 cout << "row " << i << " pivot in row " << k
488 << " colum " << j << endl;
489 }
490
491
492
493
494
495 // pivot element found in row k, check if we need to swap rows:
496 if (k != i) {
497 for (jj = 0; jj < n; jj++) {
498 Algo.int_swap(A[i * n + jj], A[k * n + jj]);
499 }
500 }
501
502
503 // now, pivot is in row i, column j :
504
505 pivot = A[i * n + j];
506 if (f_vv) {
507 cout << "pivot=" << pivot << endl;
508 }
509 pivot_inv = F->inverse(pivot);
510 if (f_vv) {
511 cout << "pivot=" << pivot << " pivot_inv=" << pivot_inv << endl;
512 }
513 if (!f_special) {
514 // make pivot to 1:
515 for (jj = 0; jj < n; jj++) {
516 A[i * n + jj] = F->mult(A[i * n + jj], pivot_inv);
517 }
518 if (f_vv) {
519 cout << "pivot=" << pivot << " pivot_inv=" << pivot_inv
520 << " made to one: " << A[i * n + j] << endl;
521 }
522 if (f_vvv) {
523 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
524 }
525 }
526
527 // do the gaussian elimination:
528
529 if (f_vv) {
530 cout << "doing elimination in column " << j << " from row "
531 << i + 1 << " down to row " << m - 1 << ":" << endl;
532 }
533 for (k = i + 1; k < m; k++) {
534 if (f_vv) {
535 cout << "k=" << k << endl;
536 }
537 z = A[k * n + j];
538 if (z == 0) {
539 continue;
540 }
541 if (f_special) {
542 f = F->mult(z, pivot_inv);
543 }
544 else {
545 f = z;
546 }
547 f = F->negate(f);
548 //A[k * n + j] = 0;
549 if (f_vv) {
550 cout << "eliminating row " << k << endl;
551 }
552 for (jj = 0; jj < n; jj++) {
553 a = A[i * n + jj];
554 b = A[k * n + jj];
555 // c := b + f * a
556 // = b - z * a if !f_special
557 // b - z * pivot_inv * a if f_special
558 c = F->mult(f, a);
559 c = F->add(c, b);
560 A[k * n + jj] = c;
561 if (f_vv) {
562 cout << A[k * n + jj] << " ";
563 }
564 }
565 if (f_vv) {
566 cout << endl;
567 }
568 if (f_vvv) {
569 cout << "A=" << endl;
570 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
571 }
572 }
573 if (f_vv) {
574 cout << "A=" << endl;
575 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
576 //print_integer_matrix(cout, A, m, n);
577 }
578 } // next j
579 if (f_complete) {
580 for (i = nb_pivots - 1; i >= 0; i--) {
581 j = pivots[i];
582 if (!f_special) {
583 a = A[i * n + j];
584 }
585 else {
586 pivot = A[i * n + j];
587 pivot_inv = F->inverse(pivot);
588 }
589 // do the gaussian elimination in the upper part:
590 for (k = i - 1; k >= 0; k--) {
591 z = A[k * n + j];
592 if (z == 0) {
593 continue;
594 }
595 //A[k * n + j] = 0;
596 for (jj = 0; jj < n; jj++) {
597 a = A[i * n + jj];
598 b = A[k * n + jj];
599 if (f_special) {
600 a = F->mult(a, pivot_inv);
601 }
602 c = F->mult(z, a);
603 c = F->negate(c);
604 c = F->add(c, b);
605 A[k * n + jj] = c;
606 }
607 } // next k
608 } // next i
609 }
610 if (f_vv) {
611 cout << endl;
612 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
613 //print_integer_matrix(cout, A, rank, n);
614 }
615 if (f_v) {
616 cout << "linear_algebra::Gauss_int_with_given_pivots done" << endl;
617 }
618 return TRUE;
619}
620
621
622
623int linear_algebra::RREF_search_pivot(int *A, int m, int n,
624 int &i, int &j, int *base_cols, int verbose_level)
625// A is a m x n matrix,
626{
627 int f_v = (verbose_level >= 1);
628 int f_vv = (verbose_level >= 2);
629 int k, jj;
631
632 if (f_v) {
633 cout << "linear_algebra::RREF_search_pivot" << endl;
634 }
635 if (f_vv) {
636 cout << "linear_algebra::RREF_search_pivot matrix:" << endl;
637 Int_vec_print_integer_matrix_width(cout, A, m, n, n, 5);
638 //print_tables();
639 }
640 for (; j < n; j++) {
641 if (f_vv) {
642 cout << "j=" << j << endl;
643 }
644 /* search for pivot element: */
645 for (k = i; k < m; k++) {
646 if (A[k * n + j]) {
647 if (f_vv) {
648 cout << "i=" << i << " pivot found in "
649 << k << "," << j << endl;
650 }
651 // pivot element found:
652 if (k != i) {
653 for (jj = j; jj < n; jj++) {
654 Algo.int_swap(A[i * n + jj], A[k * n + jj]);
655 }
656 }
657 break;
658 } // if != 0
659 } // next k
660
661 if (k == m) { // no pivot found
662 if (f_vv) {
663 cout << "no pivot found" << endl;
664 }
665 continue; // increase j, leave i constant
666 }
667
668 if (f_vv) {
669 cout << "row " << i << " pivot in row "
670 << k << " colum " << j << endl;
671 }
672 base_cols[i] = j;
673 return TRUE;
674 } // next j
675 return FALSE;
676}
677
678void linear_algebra::RREF_make_pivot_one(int *A, int m, int n,
679 int &i, int &j, int *base_cols, int verbose_level)
680// A is a m x n matrix,
681{
682 int f_v = (verbose_level >= 1);
683 int f_vv = (verbose_level >= 2);
684 int pivot, pivot_inv;
685 int jj;
686
687 if (f_v) {
688 cout << "linear_algebra::RREF_make_pivot_one" << endl;
689 }
690 pivot = A[i * n + j];
691 if (f_vv) {
692 cout << "pivot=" << pivot << endl;
693 }
694 //pivot_inv = inv_table[pivot];
695 pivot_inv = F->inverse(pivot);
696 if (f_vv) {
697 cout << "pivot=" << pivot << " pivot_inv="
698 << pivot_inv << endl;
699 }
700 // make pivot to 1:
701 for (jj = j; jj < n; jj++) {
702 A[i * n + jj] = F->mult(A[i * n + jj], pivot_inv);
703 }
704 if (f_vv) {
705 cout << "pivot=" << pivot << " pivot_inv=" << pivot_inv
706 << " made to one: " << A[i * n + j] << endl;
707 }
708 if (f_v) {
709 cout << "linear_algebra::RREF_make_pivot_one done" << endl;
710 }
711}
712
713
714void linear_algebra::RREF_elimination_below(int *A, int m, int n,
715 int &i, int &j, int *base_cols, int verbose_level)
716// A is a m x n matrix,
717{
718 int f_v = (verbose_level >= 1);
719 int f_vv = (verbose_level >= 2);
720 int k, jj, z, f, a, b, c;
721
722 if (f_v) {
723 cout << "linear_algebra::RREF_elimination_below" << endl;
724 }
725 for (k = i + 1; k < m; k++) {
726 if (f_vv) {
727 cout << "k=" << k << endl;
728 }
729 z = A[k * n + j];
730 if (z == 0) {
731 continue;
732 }
733 f = z;
734 f = F->negate(f);
735 A[k * n + j] = 0;
736 if (f_vv) {
737 cout << "eliminating row " << k << endl;
738 }
739 for (jj = j + 1; jj < n; jj++) {
740 a = A[i * n + jj];
741 b = A[k * n + jj];
742 // c := b + f * a
743 // = b - z * a if !f_special
744 // b - z * pivot_inv * a if f_special
745 c = F->mult(f, a);
746 c = F->add(c, b);
747 A[k * n + jj] = c;
748 if (f_vv) {
749 cout << A[k * n + jj] << " ";
750 }
751 }
752 }
753 i++;
754 if (f_v) {
755 cout << "linear_algebra::RREF_elimination_below done" << endl;
756 }
757}
758
759void linear_algebra::RREF_elimination_above(int *A, int m, int n,
760 int i, int *base_cols, int verbose_level)
761// A is a m x n matrix,
762{
763 int f_v = (verbose_level >= 1);
764 int j, k, jj, z, a, b, c;
765
766 if (f_v) {
767 cout << "linear_algebra::RREF_elimination_above" << endl;
768 }
769 j = base_cols[i];
770 a = A[i * n + j];
771 // do the gaussian elimination in the upper part:
772 for (k = i - 1; k >= 0; k--) {
773 z = A[k * n + j];
774 if (z == 0) {
775 continue;
776 }
777 A[k * n + j] = 0;
778 for (jj = j + 1; jj < n; jj++) {
779 a = A[i * n + jj];
780 b = A[k * n + jj];
781 c = F->mult(z, a);
782 c = F->negate(c);
783 c = F->add(c, b);
784 A[k * n + jj] = c;
785 }
786 } // next k
787 if (f_v) {
788 cout << "linear_algebra::RREF_elimination_above done" << endl;
789 }
790}
791
792
793}}}
794
#define Int_vec_print_integer_matrix(A, B, C, D)
Definition: foundations.h:690
#define Int_vec_print_integer_matrix_width(A, B, C, D, E, F)
Definition: foundations.h:691
#define TRUE
Definition: foundations.h:231
#define FALSE
Definition: foundations.h:234
#define Int_vec_print(A, B, C)
Definition: foundations.h:685
int Gauss_int(int *A, int f_special, int f_complete, int *base_cols, int f_P, int *P, int m, int n, int Pn, int q, int *add_table, int *mult_table, int *negate_table, int *inv_table, int verbose_level)
Definition: global.cpp:1716
the orbiter library for the classification of combinatorial objects