1 //==============================================================================================
2 //
3 // This file is part of LiDIA --- a library for computational number theory
4 //
5 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved.
6 //
7 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
8 //
9 //----------------------------------------------------------------------------------------------
10 //
11 // $Id$
12 //
13 // Author : Michael Jacobson (MJ)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18
19 #ifdef HAVE_CONFIG_H
20 # include "config.h"
21 #endif
22 #include "LiDIA/bigint_matrix.h"
23 #include <cstdlib>
24 #include <cstdio>
25 #include <fstream>
26 #include <unistd.h>
27
28
29
30 #ifdef LIDIA_NAMESPACE
31 namespace LiDIA {
32 #endif
33
34
35
trans_matrix()36 trans_matrix::trans_matrix()
37 {
38 static int curr = 0;
39
40 files = false;
41 size = 0;
42 TR.set_mode(EXPAND);
43 index = curr;
44 last_cols = 0;
45 notone_idx = -1;
46 notone_rows = 0;
47 diag_element = 1;
48
49 adj_files = false;
50 adj_index = -1;
51
52 DET_div.assign_zero();
53 det_idx = -1;
54
55 ++curr;
56 }
57
58
~trans_matrix()59 trans_matrix::~trans_matrix()
60 {
61 kill();
62 }
63
64
65 //
66 // trans_matrix::kill()
67 //
68 // Deletes any temporary files used by the trans_matrix, and reinitializes
69 // all members.
70 //
71
72 void
kill()73 trans_matrix::kill()
74 {
75 if (files) {
76 lidia_size_t i;
77 char fname[50];
78
79 for (i = 0; i < size; ++i) {
80 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, i);
81 std::remove(fname);
82 }
83 }
84
85 if (adj_files)
86 ADJ.kill();
87
88 size = 0;
89 TR.set_size(0);
90 last_cols = 0;
91 notone_idx = -1;
92 notone_rows = 0;
93 diag_element = 1;
94
95 adj_index = -1;
96 }
97
98
99
100 //
101 // trans_matrix::set_mode(int)
102 //
103 // Sets the storage mode of the trans_matrix:
104 // 0 - no files
105 // 1 - components stored in files
106 // 2 - components stored in files, adjoint also stored in files
107 //
108
109 void
set_mode(int mode)110 trans_matrix::set_mode(int mode)
111 {
112 files = adj_files = false;
113
114 if (mode > 0)
115 files = true;
116 if (mode > 1)
117 adj_files = true;
118 }
119
120
121
122 //
123 // trans_matrix::get_adjoint_mode()
124 //
125 // Returns true if the adjoint is stored in its own set of files, false
126 // otherwise.
127 //
128
129 bool
get_adjoint_mode()130 trans_matrix::get_adjoint_mode()
131 {
132 return adj_files;
133 }
134
135
136
137 //
138 // trans_matrix::touch_files
139 //
140 // Touches all the temporary files of the trans_matrix in order to give
141 // them current date stamps (i.e., so they aren't deleted by the system
142 // before we're finished).
143 //
144
145 void
touch_files()146 trans_matrix::touch_files()
147 {
148 // touch all temporary files in order to give current date stamps
149 if (files) {
150 register lidia_size_t i;
151 char command[100];
152
153 for (i = 0; i < size; ++i) {
154 if (i == adj_index)
155 ADJ.touch_files();
156 else {
157 sprintf(command, "touch /tmp/%ldQO_TRANS_%i.%i", (long)getpid(), index, i);
158 std::system(command);
159 }
160 }
161 }
162 }
163
164
165
166 //
167 // trans_matrix::remove_columns(lidia_dize_t *)
168 //
169 // Removes the selected columns from the trans_matrix.
170 //
171
172 void
remove_columns(lidia_size_t * rcols)173 trans_matrix::remove_columns(lidia_size_t *rcols)
174 {
175 if (files) {
176 std::ofstream out;
177 std::ifstream in;
178 char fname[50];
179 matrix< bigint > tran;
180
181 tran.set_storage_mode(matrix_flags::sparse_representation);
182 tran.set_orientation(matrix_flags::row_oriented);
183 tran.set_print_mode(LIDIA_MODE);
184
185 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, size-1);
186 in.open(fname);
187 if (in.fail()) {
188 std::cout << "ERROR (remove_columns)! ";
189 std::cout << "Can't open transformation matrix file for input!!!" << std::endl;
190 }
191 in >> tran;
192 in.close();
193
194 tran.remove_columns(rcols);
195
196 out.open(fname);
197 if (out.fail()) {
198 std::cout << "ERROR (remove)! ";
199 std::cout << "Can't open transformation matrix file for output!!!" << std::endl;
200 }
201 out << tran;
202 out.close();
203
204 last_cols = tran.columns;
205 }
206 else {
207 TR[size-1].remove_columns(rcols);
208 last_cols = TR[size-1].columns;
209 }
210 }
211
212
213
214 //
215 // trans_matrix::set_no_of_columns(lidia_size_t)
216 //
217 // Sets the number of columns of the trans_matrix, either by adding or
218 // deleting.
219 //
220
221 void
set_no_of_columns(lidia_size_t ncols)222 trans_matrix::set_no_of_columns(lidia_size_t ncols)
223 {
224 lidia_size_t i, j, k, ii, old_rows;
225 bigint newel;
226
227 if (size == 0) {
228 cols = ncols;
229 return;
230 }
231
232 if (files) {
233 std::ofstream out;
234 std::ifstream in;
235 char fname[50];
236 matrix< bigint > tran;
237
238 tran.set_storage_mode(matrix_flags::sparse_representation);
239 tran.set_orientation(matrix_flags::row_oriented);
240 tran.set_print_mode(LIDIA_MODE);
241
242 if (ncols > cols) {
243 k = ncols - cols;
244 for (i = 0; i < size; ++i) {
245 if (i == adj_index) {
246 ADJ.resize(ncols, ncols);
247 last_cols = ncols;
248 }
249 else {
250 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, i);
251 in.open(fname);
252 if (in.fail()) {
253 std::cout << "ERROR (set_no_of_columns)! ";
254 std::cout << "Can't open trans matrix file for input (1)!!!" << std::endl;
255 }
256 in >> tran;
257 in.close();
258
259 old_rows = tran.get_no_of_rows();
260 ii = tran.get_no_of_columns();
261 if (i == notone_idx)
262 newel = diag_element;
263 else
264 newel = 1;
265 tran.resize(old_rows+k, ii+k);
266 for (j = old_rows; j < old_rows+k; ++j) {
267 tran.sto(j, ii, newel);
268 ++ii;
269 }
270
271 out.open(fname);
272 if (out.fail()) {
273 std::cout << "ERROR (set_no_of_columns)! ";
274 std::cout << "Can't open trans matrix file for output (1)!!!" << std::endl;
275 }
276 out << tran;
277 out.close();
278
279 last_cols = tran.columns;
280 }
281 }
282 }
283 else {
284 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, size-1);
285 in.open(fname);
286 if (in.fail()) {
287 std::cout << "ERROR (set_no_of_columns)! ";
288 std::cout << "Can't open trans matrix file for input (2)!!!" << std::endl;
289 }
290 in >> tran;
291 in.close();
292
293 tran.set_no_of_columns(ncols);
294
295 out.open(fname);
296 if (out.fail()) {
297 std::cout << "ERROR (set_no_of_columns)! ";
298 std::cout << "Can't open trans matrix file for output (2)!!!" << std::endl;
299 }
300 out << tran;
301 out.close();
302
303 last_cols = tran.columns;
304 }
305 }
306 else {
307 if (ncols > cols) {
308 k = ncols - cols;
309 for (i = 0; i < size; ++i) {
310 old_rows = TR[i].get_no_of_rows();
311 ii = TR[i].get_no_of_columns();
312 if (i == notone_idx)
313 newel = diag_element;
314 else
315 newel = 1;
316 TR[i].resize(old_rows+k, ii+k);
317 for (j = old_rows; j < old_rows+k; ++j) {
318 TR[i].sto(j, ii, newel);
319 ++ii;
320 }
321 }
322 }
323 else
324 TR[size-1].set_no_of_columns(ncols);
325
326 last_cols = TR[size-1].columns;
327 }
328
329 cols = ncols;
330 }
331
332
333
334 //
335 // trans_matrix::get_no_of_columns()
336 //
337 // Returns the number of columns
338 //
339
340 lidia_size_t
get_no_of_columns()341 trans_matrix::get_no_of_columns()
342 {
343 return cols;
344 }
345
346
347
348 //
349 // trans_matrix::get_size()
350 //
351 // Returns the numberer of component matrices
352 //
353
354 lidia_size_t
get_size()355 trans_matrix::get_size()
356 {
357 return size;
358 }
359
360
361
362 //
363 // trans_matrix::set_size()
364 //
365 // Sets the numberer of component matrices
366 //
367
368 void
set_size(lidia_size_t nsize)369 trans_matrix::set_size(lidia_size_t nsize)
370 {
371 if (files) {
372 lidia_size_t i;
373 std::ofstream out;
374 char fname[50];
375
376 if (nsize < size) {
377 for (i = nsize; i < size; ++i) {
378 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, i);
379 std::remove(fname);
380 }
381 }
382 else if (nsize > size) {
383 for (i = size; i < nsize; ++i) {
384 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, i);
385 out.open(fname);
386 out.close();
387 }
388 }
389 }
390 else
391 TR.set_size(nsize);
392
393 size = nsize;
394 }
395
396
397
398 //
399 // trans_matrix::get_matrix(matrix <bigint> &, lidia_size_t)
400 //
401 // Loads (from memory or disk) the component matrix with index idx into
402 // the matrix tran
403 //
404
405 void
get_matrix(matrix<bigint> & tran,lidia_size_t idx)406 trans_matrix::get_matrix(matrix< bigint > & tran, lidia_size_t idx)
407 {
408 if (idx >= size) {
409 std::cout << "Only " << size << "matrices (input " << idx << ")" << std::endl;
410 return;
411 }
412
413 unsigned long osmode = tran.get_storage_mode();
414 unsigned long oorient = tran.get_orientation();
415 unsigned long opmode = tran.get_print_mode();
416
417 if (files) {
418 if (idx == adj_index) {
419 std::cout << "Can't get adjoint from file (not implemented (or needed)" << std::endl;
420 return;
421 }
422
423 std::ifstream in;
424 char fname[50];
425
426 tran.set_storage_mode(matrix_flags::sparse_representation);
427 tran.set_orientation(matrix_flags::row_oriented);
428 tran.set_print_mode(LIDIA_MODE);
429
430 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, idx);
431 in.open(fname);
432 if (in.fail()) {
433 std::cout << "ERROR (get_matrix)! ";
434 std::cout << "Can't open transformation matrix file for input!!!" << std::endl;
435 }
436 in >> tran;
437 in.close();
438
439 tran.set_storage_mode(osmode);
440 tran.set_orientation(oorient);
441 tran.set_print_mode(opmode);
442 }
443 else
444 tran.assign(TR[idx]);
445 }
446
447
448
449 //
450 // trans_matrix::store_matrix(matrix <bigint> &)
451 //
452 // Stores the matrix tran as the next component matrix in the trans_matrix,
453 // either in memory or to disk as indicated by the mode.
454 //
455
456 void
store_matrix(matrix<bigint> & tran)457 trans_matrix::store_matrix(matrix< bigint > & tran)
458 {
459 std::ofstream out;
460 lidia_size_t i, TRcols = tran.columns, TRrows = tran.rows;
461 char fname[50];
462 unsigned long osmode = tran.get_storage_mode();
463 unsigned long oorient = tran.get_orientation();
464 unsigned long opmode = tran.get_print_mode();
465
466 tran.set_storage_mode(matrix_flags::sparse_representation);
467 tran.set_orientation(matrix_flags::row_oriented);
468 if (size == 0) {
469 tran.resize(cols, cols);
470 for (i = TRcols; i < cols; ++i)
471 tran.sto(i, i, bigint(1));
472 last_cols = cols;
473 }
474 else {
475 tran.resize(last_cols, last_cols);
476 for (i = TRcols; i < last_cols; ++i)
477 tran.sto(i, i, bigint(1));
478 }
479
480 if (files) {
481 tran.set_print_mode(LIDIA_MODE);
482 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, size);
483 out.open(fname);
484 if (out.fail()) {
485 std::cout << "ERROR (store_matrix)! ";
486 std::cout << "Can't open transformation matrix file for output!!!" << std::endl;
487 }
488 out << tran;
489 out.close();
490 }
491 else
492 TR[size].assign(tran);
493
494 ++size;
495
496 tran.reset();
497 tran.set_storage_mode(osmode);
498 tran.set_orientation(oorient);
499 tran.set_print_mode(opmode);
500 tran.resize(TRrows, TRcols);
501
502 touch_files();
503 }
504
505
506
507 //
508 // trans_matrix::store_matrix(matrix <bigint> &, bigint &)
509 //
510 // Stores the matrix tran as the next component matrix in the trans_matrix,
511 // either in memory or to disk as indicated by the mode. DET is the multiple
512 // by which the additional diagonal elements must be multiplied in the case
513 // that AT = D HNF(A).
514 //
515
516 void
store_matrix(matrix<bigint> & tran,bigint & DET)517 trans_matrix::store_matrix(matrix< bigint > & tran, bigint & DET)
518 {
519 std::ofstream out;
520 lidia_size_t i, TRcols = tran.columns, TRrows = tran.rows;
521 char fname[50];
522 unsigned long osmode = tran.get_storage_mode();
523 unsigned long oorient = tran.get_orientation();
524 unsigned long opmode = tran.get_print_mode();
525
526 notone_idx = size;
527 diag_element.assign(DET);
528
529 tran.set_storage_mode(matrix_flags::sparse_representation);
530 tran.set_orientation(matrix_flags::row_oriented);
531 if (size == 0) {
532 tran.resize(cols, cols);
533 for (i = TRcols; i < cols; ++i)
534 tran.sto(i, i, DET);
535 last_cols = cols;
536 }
537 else {
538 tran.resize(last_cols, last_cols);
539 for (i = TRcols; i < last_cols; ++i)
540 tran.sto(i, i, DET);
541 }
542
543
544 if (files) {
545 tran.set_print_mode(LIDIA_MODE);
546 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, size);
547 out.open(fname);
548 if (out.fail()) {
549 std::cout << "ERROR (store_matrix_DET)! ";
550 std::cout << "Can't open transformation matrix file!!! for output" << std::endl;
551 }
552 out << tran;
553 out.close();
554 }
555 else
556 TR[size].assign(tran);
557
558 ++size;
559
560 tran.reset();
561 tran.set_storage_mode(osmode);
562 tran.set_orientation(oorient);
563 tran.set_print_mode(opmode);
564 tran.resize(TRrows, TRcols);
565
566 touch_files();
567 }
568
569
570
571 //
572 // trans_matrix::store_adjoint
573 //
574 // Sets the component index of the adjoint
575 //
576
577 void
store_adjoint()578 trans_matrix::store_adjoint()
579 {
580 ADJ.resize(cols, cols);
581 adj_index = size;
582 ++size;
583 }
584
585
586
587 //
588 // trans_matrix::store_det
589 //
590 // Sets the component index of the matrix whose diagonal elements must be
591 // multiplied by DET
592 //
593
594 void
store_det(bigint & DET)595 trans_matrix::store_det(bigint & DET)
596 {
597 DET_div.assign(DET);
598 det_idx = size;
599 }
600
601
602
603 //
604 // trans_matrix::get_column_vector
605 //
606 // Computes column idx of the trans_matrix
607 //
608
609 void
get_column_vector(math_vector<bigint> & vec,lidia_size_t idx)610 trans_matrix::get_column_vector(math_vector< bigint > & vec, lidia_size_t idx)
611 {
612 debug_handler("trans_matrix", "make_column");
613
614 lidia_size_t i, j, k, last;
615 math_vector< bigint > tmp;
616 bigint ajk;
617
618 tmp.set_mode(EXPAND);
619
620 if (files) {
621 std::ifstream in;
622 char fname[50];
623 matrix< bigint > tran;
624
625 tran.set_storage_mode(matrix_flags::sparse_representation);
626 tran.set_orientation(matrix_flags::row_oriented);
627 tran.set_print_mode(LIDIA_MODE);
628
629 // get last matrix and the desired column
630 last = size-1;
631
632 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, last);
633 in.open(fname);
634 if (in.fail()) {
635 std::cout << "ERROR (get_column_vector)! ";
636 std::cout << "Can't open trans matrix file for input (1)!!!" << std::endl;
637 }
638 in >> tran;
639 in.close();
640
641 tran.get_column_vector(vec, idx);
642 tran.reset();
643
644 // evaluate product
645 for (i = last-1; i >= 0; --i) {
646 if (i == adj_index) {
647 ADJ.multiply(tmp, vec);
648 vec = tmp;
649 }
650 else {
651 // get matrix from file
652 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, i);
653 in.open(fname);
654 if (in.fail()) {
655 std::cout << "ERROR (get_column_vector)! ";
656 std::cout << "Can't open trans matrix file for input (2)!!!" << std::endl;
657 }
658 in >> tran;
659 in.close();
660
661 // perform non-trivial multiplication (multiply(vec,TR[i],vec))
662 for (j = 0; j < tran.rows; ++j) {
663 tmp[j] = 0;
664 for (k = 0; k < tran.value_counter[j]; ++k) {
665 ajk.assign(tran.value[j][k]);
666 LiDIA::multiply(ajk, ajk, vec[tran.index[j][k]]);
667 add(tmp[j], tmp[j], ajk);
668 }
669 }
670
671 vec = tmp;
672 tran.reset();
673 }
674
675 if (i == det_idx)
676 LiDIA::divide(vec, vec, DET_div);
677 }
678 }
679 else {
680 last = size-1;
681 TR[last].get_column_vector(vec, idx);
682
683 for (i = last-1; i >= 0; --i) {
684 // perform non-trivial multiplication (multiply(vec,TR[i],vec))
685 for (j = 0; j < TR[i].rows; ++j) {
686 tmp[j] = 0;
687 for (k = 0; k < TR[i].value_counter[j]; ++k) {
688 ajk.assign(TR[i].value[j][k]);
689 LiDIA::multiply(ajk, ajk, vec[TR[i].index[j][k]]);
690 add(tmp[j], tmp[j], ajk);
691 }
692 }
693
694 vec = tmp;
695 if (i == det_idx)
696 LiDIA::divide(vec, vec, DET_div);
697 }
698 }
699 }
700
701
702
703 //
704 // trans_matrix::get_submatrix
705 //
706 // Computes the submatrix containing the given columns of the trans_matrix
707 //
708
709 void
get_submatrix(matrix<bigint> & V,lidia_size_t * idx)710 trans_matrix::get_submatrix(matrix< bigint > & V, lidia_size_t *idx)
711 {
712 debug_handler("trans_matrix", "get_submatrix");
713
714 if (!idx)
715 return;
716
717 lidia_size_t i, j, k, l, last;
718 math_vector< bigint > vec, tmp;
719 bigint ajk, val;
720
721 vec.set_mode(EXPAND);
722 tmp.set_mode(EXPAND);
723
724 if (files) {
725 std::ifstream in;
726 char fname[50];
727 matrix< bigint > tran;
728
729 tran.set_storage_mode(matrix_flags::sparse_representation);
730 tran.set_orientation(matrix_flags::row_oriented);
731 tran.set_print_mode(LIDIA_MODE);
732
733 // get last matrix and the desired column
734 last = size-1;
735
736 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, last);
737 in.open(fname);
738 if (in.fail()) {
739 std::cout << "ERROR (get_submatrix)! ";
740 std::cout << "Can't open trans matrix file for input (1)!!!" << std::endl;
741 }
742 in >> tran;
743 in.close();
744
745 V.resize(tran.rows, idx[0]);
746 for (i = 0; i < idx[0]; ++i) {
747 tran.get_column_vector(vec, idx[i+1]);
748 V.sto_column_vector(vec, tran.rows, i);
749 }
750 tran.reset();
751
752 // evaluate product
753 for (i = last-1; i >= 0; --i) {
754 if (i == adj_index) {
755 matrix< bigint > B;
756 ADJ.multiply(B, V);
757 V = B;
758 }
759 else {
760 // get matrix from file
761 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, i);
762 in.open(fname);
763 if (in.fail()) {
764 std::cout << "ERROR (get_submatrix)! ";
765 std::cout << "Can't open trans matrix file for input (2)!!!" << std::endl;
766 }
767 in >> tran;
768 in.close();
769
770 // perform non-trivial multiplication (multiply(vec,TR[i],vec))
771 for (l = 0; l < V.columns; ++l) {
772 V.get_column_vector(vec, l);
773
774 for (j = 0; j < tran.rows; ++j) {
775 tmp[j] = 0;
776 for (k = 0; k < tran.value_counter[j]; ++k) {
777 ajk.assign(tran.value[j][k]);
778 LiDIA::multiply(ajk, ajk, vec[tran.index[j][k]]);
779 add(tmp[j], tmp[j], ajk);
780 }
781 }
782
783 V.sto_column_vector(tmp, V.rows, l);
784 }
785
786 tran.reset();
787 }
788
789 if (i == det_idx)
790 LiDIA::divide(V, V, DET_div);
791 }
792 }
793 else {
794 last = size-1;
795 V.resize(TR[last].rows, idx[0]);
796 for (i = 0; i < idx[0]; ++i) {
797 TR[last].get_column_vector(vec, idx[i+1]);
798 V.sto_column_vector(vec, TR[last].rows, i);
799 }
800
801 for (i = last-1; i >= 0; --i) {
802 // perform non-trivial multiplication (multiply(V,TR[i],V))
803 for (l = 0; l < V.columns; ++l) {
804 V.get_column_vector(vec, l);
805
806 for (j = 0; j < TR[i].rows; ++j) {
807 tmp[j] = 0;
808 for (k = 0; k < TR[i].value_counter[j]; ++k) {
809 ajk.assign(TR[i].value[j][k]);
810 LiDIA::multiply(ajk, ajk, vec[TR[i].index[j][k]]);
811 add(tmp[j], tmp[j], ajk);
812 }
813 }
814
815 V.sto_column_vector(tmp, V.rows, l);
816 }
817
818 if (i == det_idx)
819 LiDIA::divide(V, V, DET_div);
820 }
821 }
822 }
823
824
825
826 //
827 // trans_matrix::multiply()
828 //
829 // Computes the product of the trans_matrix with the vector x, i.e., Tx.
830 //
831
832 math_vector< bigint >
multiply(math_vector<bigint> & x)833 trans_matrix::multiply(math_vector< bigint > & x)
834 {
835 debug_handler("trans_matrix", "multiply");
836
837 lidia_size_t i, j, k, last;
838 math_vector< bigint > tmp, vec;
839 bigint ajk;
840
841 vec.set_mode(EXPAND);
842 tmp.set_mode(EXPAND);
843 vec = x;
844
845 if (files) {
846 std::ifstream in;
847 char fname[50];
848 matrix< bigint > tran;
849
850 tran.set_storage_mode(matrix_flags::sparse_representation);
851 tran.set_orientation(matrix_flags::row_oriented);
852 tran.set_print_mode(LIDIA_MODE);
853
854 // get last matrix
855 last = size;
856
857 // evaluate product
858 for (i = last-1; i >= 0; --i) {
859 if (i == adj_index) {
860 ADJ.multiply(tmp, vec);
861 vec = tmp;
862 }
863 else {
864 // get matrix from file
865 sprintf(fname, "/tmp/%ldQO_TRANS_%d.%d", (long) getpid(), index, i);
866 in.open(fname);
867 if (in.fail()) {
868 std::cout << "ERROR (multiply)! ";
869 std::cout << "Can't open transformation matrix file for input!!!" << std::endl;
870 }
871 in >> tran;
872 in.close();
873
874 // perform non-trivial multiplication (multiply(vec,TR[i],vec))
875 for (j = 0; j < tran.rows; ++j) {
876 tmp[j] = 0;
877 for (k = 0; k < tran.value_counter[j]; ++k) {
878 ajk.assign(tran.value[j][k]);
879 LiDIA::multiply(ajk, ajk, vec[tran.index[j][k]]);
880 add(tmp[j], tmp[j], ajk);
881 }
882 }
883
884 vec = tmp;
885 tran.reset();
886 }
887
888 if (i == det_idx)
889 LiDIA::divide(vec, vec, DET_div);
890 }
891 }
892 else {
893 last = size;
894
895 for (i = last-1; i >= 0; --i) {
896 // perform non-trivial multiplication (multiply(vec,TR[i],vec))
897 for (j = 0; j < TR[i].rows; ++j) {
898 tmp[j] = 0;
899 for (k = 0; k < TR[i].value_counter[j]; ++k) {
900 ajk.assign(TR[i].value[j][k]);
901 LiDIA::multiply(ajk, ajk, vec[TR[i].index[j][k]]);
902 add(tmp[j], tmp[j], ajk);
903 }
904 }
905
906 vec = tmp;
907
908 if (i == det_idx)
909 LiDIA::divide(vec, vec, DET_div);
910 }
911 }
912
913 return vec;
914 }
915
916
917
918 #ifdef LIDIA_NAMESPACE
919 } // end of namespace LiDIA
920 #endif
921