source: git/IntegerProgramming/matrix.cc @ 3e0c94f

spielwiese
Last change on this file since 3e0c94f was 4dcfb4f, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: const Bigint&)0 git-svn-id: file:///usr/local/Singular/svn/trunk@7681 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 17.2 KB
Line 
1// matrix.cc
2
3// implementation of class matrix
4
5#ifndef MATRIX_CC
6#define MATRIX_CC
7
8#include "matrix.h"
9
10////////////// constructors and destructor //////////////////////////////////
11typedef Integer* IntegerP;
12typedef BigInt* BigIntP;
13
14matrix::matrix(const short& row_number, const short& column_number)
15    :rows(row_number),columns(column_number)
16{
17  _kernel_dimension=-2;
18  // LLL-algorithm not yet performed
19
20  // argument check
21  if((rows<=0)||(columns<=0))
22    // bad input, set "error flag"
23  {
24    cerr<<"\nWARNING: matrix::matrix(const short&, const short&):\n"
25      "argument out of range"<<endl;
26    columns=-1;
27    return;
28  }
29
30  // memory allocation and initialization
31
32  coefficients=new IntegerP[rows];
33  for(short i=0;i<rows;i++)
34    coefficients[i]=new Integer[columns];
35  for(short i=0;i<rows;i++)
36    for(short j=0;j<columns;j++)
37      coefficients[i][j]=0;
38}
39
40matrix::matrix(const short& row_number, const short& column_number,
41               Integer** entries)
42    :rows(row_number),columns(column_number)
43{
44  _kernel_dimension=-2;
45  // LLL-algorithm not yet performed
46
47  // argument check
48  if((rows<=0)||(columns<=0))
49    // bad input, set "error flag"
50  {
51    cerr<<"\nWARNING: matrix::matrix(const short&, const short&, Integr**):\n"
52      "argument out of range"<<endl;
53    columns=-1;
54    return;
55  }
56
57  // memory allocation and initialization
58
59  coefficients=new IntegerP[rows];
60  for(short i=0;i<rows;i++)
61    coefficients[i]=new Integer[columns];
62  for(short i=0;i<rows;i++)
63    for(short j=0;j<columns;j++)
64      coefficients[i][j]=entries[i][j];
65  // coefficients[i] is the i-th row vector
66}
67
68
69
70
71matrix::matrix(ifstream& input)
72{
73  _kernel_dimension=-2;
74  // LLL-algorithm not yet performed
75
76  input>>rows;
77  if(!input)
78    // input failure, set "error flag"
79  {
80    cerr<<"\nWARNING: matrix::matrix(ifstream&): input failure"<<endl;
81    columns=-2;
82    return;
83  }
84
85  input>>columns;
86  if(!input)
87    // input failure, set "error flag"
88  {
89    cerr<<"\nWARNING: matrix::matrix(ifstream&): input failure"<<endl;
90    columns=-2;
91    return;
92  }
93
94  if((rows<=0)||(columns<=0))
95    // bad input, set "error flag"
96  {
97    cerr<<"\nWARNING: matrix::matrix(ifstream&): bad input"<<endl;
98    columns=-1;
99    return;
100  }
101
102  coefficients=new IntegerP[rows];
103  for(short i=0;i<rows;i++)
104    coefficients[i]=new Integer[columns];
105  for(short i=0;i<rows;i++)
106    for(short j=0;j<columns;j++)
107    {
108      input>>coefficients[i][j];
109      if(!input)
110        // bad input, set error flag
111      {
112        cerr<<"\nWARNING: matrix::matrix(ifstream&): input failure"<<endl;
113        columns=-2;
114        return;
115      }
116    }
117}
118
119
120
121
122matrix::matrix(const short& m, const short& n, ifstream& input)
123{
124  _kernel_dimension=-2;
125  // LLL-algorithm not yet performed
126
127  // argument check
128  if((m<=0) || (n<=0))
129    // bad input, set "error flag"
130  {
131    cerr<<"\nWARNING: matrix::matrix(const short&, const short&, ifstream&):\n"
132      "argument out of range"<<endl;
133    columns=-1;
134    return;
135  }
136
137  rows=m;
138  columns=n;
139
140  // memory allocation and initialization
141
142  coefficients=new IntegerP[rows];
143  for(short i=0;i<rows;i++)
144    coefficients[i]=new Integer[columns];
145  for(short i=0;i<rows;i++)
146    for(short j=0;j<columns;j++)
147    {
148      input>>coefficients[i][j];
149      if(!input)
150        // bad input, set error flag
151      {
152        columns=-2;
153        return;
154      }
155    }
156}
157
158
159
160
161matrix::matrix(const matrix& A)
162    :rows(A.rows),columns(A.columns),_kernel_dimension(A._kernel_dimension)
163{
164
165  if(columns<0)
166  {
167    cerr<<"\nWARNING: matrix::matrix(const matrix&):\n"
168      "Building a matrix from a corrupt one"<<endl;
169    return;
170  }
171
172  // memory allocation and initialization (also for H)
173
174  coefficients=new IntegerP[rows];
175  for(short i=0;i<rows;i++)
176    coefficients[i]=new Integer[columns];
177  for(short i=0;i<rows;i++)
178    for(short j=0;j<columns;j++)
179      coefficients[i][j]=A.coefficients[i][j];
180
181  if(_kernel_dimension>0)
182  {
183    H=new BigIntP[_kernel_dimension];
184    for(short k=0;k<_kernel_dimension;k++)
185      H[k]=new BigInt[columns];
186    for(short k=0;k<_kernel_dimension;k++)
187      for(short j=0;j<columns;j++)
188        H[k][j]=(A.H)[k][j];
189  }
190}
191
192
193
194
195matrix::~matrix()
196{
197  for(short i=0;i<rows;i++)
198    delete[] coefficients[i];
199  delete[] coefficients;
200
201  if(_kernel_dimension>0)
202    // LLL-algorithm performed
203  {
204    for(short i=0;i<_kernel_dimension;i++)
205      delete[] H[i];
206    delete[] H;
207  }
208}
209
210
211
212
213//////////////////// object properties //////////////////////////////////////
214
215
216
217
218BOOLEAN matrix::is_nonnegative() const
219{
220  for(short i=0;i<rows;i++)
221    for(short j=0;j<columns;j++)
222      if(coefficients[i][j]<0)
223        return FALSE;
224  return TRUE;
225}
226
227
228
229short matrix::error_status() const
230{
231  if(columns<0)
232    return columns;
233  else
234    return 0;
235}
236
237
238
239short matrix::row_number() const
240{
241  return rows;
242}
243
244
245
246short matrix::column_number() const
247{
248  return columns;
249}
250
251
252
253
254////////// special routines for the IP-algorithms /////////////////////////
255
256
257
258
259short matrix::LLL_kernel_basis()
260{
261
262  // copy the column vectors of the actual matrix
263  // (They are modified by the LLL-algorithm!)
264  BigInt** b=new BigIntP[columns];
265  for(short n=0;n<columns;n++)
266    b[n]=new BigInt[rows];
267  for(short n=0;n<columns;n++)
268    for(short m=0;m<rows;m++)
269      b[n][m]=coefficients[m][n];
270
271  // compute a LLL-reduced basis of the relations of b[0],...,b[columns-1]
272  _kernel_dimension=relations(b,columns,rows,H);
273
274  // The kernel lattice basis is now stored in the member H (vectors
275  // H[0],...,H[_kernel_dimension-1]).
276
277  // delete auxiliary vectors
278  for(short n=0;n<columns;n++)
279    delete[] b[n];
280  delete[] b;
281
282  return _kernel_dimension;
283}
284
285
286
287
288short matrix::compute_nonzero_kernel_vector()
289{
290
291  if(_kernel_dimension==-2)
292    // lattice basis not yet computed
293    LLL_kernel_basis();
294
295  if(_kernel_dimension==-1)
296  {
297    cerr<<"\nWARNING: short matrix::compute_non_zero_kernel_vector(BigInt*&):"
298      "\nerror in kernel basis, cannot compute the desired vector"<<endl;
299    return 0;
300  }
301
302  if(_kernel_dimension==0)
303  {
304    cerr<<"\nWARNING: short matrix::compute_non_zero_kernel_vector(BigInt*&): "
305      "\nkernel dimension is zero"<<endl;
306    return 0;
307  }
308
309  // Now, the kernel dimension is positive.
310
311  BigInt *M=new BigInt[_kernel_dimension];
312  // M stores a number by which the algorithm decides which vector to
313  // take next.
314
315
316// STEP 1: Determine the vector with the least zero components (if it is not
317// unique, choose the smallest).
318
319  // determine number of zero components
320  for(short i=0;i<_kernel_dimension;i++)
321  {
322    M[i]=0;
323    for(short j=0;j<columns;j++)
324      if(H[i][j]==BigInt(0))
325        M[i]++;
326  }
327
328  // determine minimal number of zero components
329  BigInt min=columns;
330  // columns is an upper bound (not reached because the kernel basis cannot
331  // contain the zero vector)
332  for(short i=0;i<_kernel_dimension;i++)
333    if(M[i]<min)
334      min=M[i];
335
336  // add the square of the norm to the vectors with the least zero components
337  // and discard the others (the norm computation is why we have chosen the
338  // M[i] to be BigInts)
339  for(short i=0;i<_kernel_dimension;i++)
340    if(M[i]!=min)
341      M[i]=-1;
342    else
343      for(short j=0;j<columns;j++)
344        M[i]+=H[i][j]*H[i][j];
345  // As the lattice basis does not contain the zero vector, at least one M[i]
346  // is positive!
347
348  // determine the start vector, i.e. the one with least zero components, but
349  // smallest possible (euclidian) norm
350  short min_index=-1;
351  for(short i=0;i<_kernel_dimension;i++)
352    if(M[i]>BigInt(0))
353      if(min_index==-1)
354        min_index=i;
355      else
356        if(M[i]<M[min_index])
357          min_index=i;
358
359  // Now, H[min_index] is the vector to be transformed into a nonnegative one.
360  // For a better overview, it is swapped with the first vector
361  // (only pointers).
362
363  if(min_index!=0)
364  {
365    BigInt* swap=H[min_index];
366    H[min_index]=H[0];
367    H[0]=swap;
368  }
369
370
371// Now construct the desired vector.
372// This is done by adding a linear combination of
373// H[1],...,H[_kernel_dimension-1] to H[0]. It is important that the final
374// result, written as a linear combination of
375// H[0],...,H[_kernel_dimension-1], has coefficient 1 or -1 at H[0]
376// (to make sure that it is together with H[1],...,H[_kernel_dimension]
377// still a  l a t t i c e   basis).
378
379  for(short current_position=1;current_position<columns;current_position++)
380    // in fact, this loop will terminate before the condition in the
381    // for-statement is satisfied...
382  {
383
384
385// STEP 2: Nonnegative vector already found?
386
387    BOOLEAN found=TRUE;
388    for(short j=0;j<columns;j++)
389      if(H[0][j]==BigInt(0))
390        found=FALSE;
391
392    if(found==TRUE)
393      // H[0] has only positive entries,
394      return 1;
395    // else there are further zero components
396
397
398// STEP 3: Can a furhter zero component be "eliminated"?
399// If this is the case, find a basis vector that can do this.
400
401    // determine number of components in each remaining vector that are zero
402    // in the vector itself as well as in the already constructed vector
403    for(short i=current_position;i<_kernel_dimension;i++)
404      M[i]=0;
405
406    short remaining_zero_components=0;
407
408    for(short j=0;j<columns;j++)
409      if(H[0][j]==BigInt(0))
410      {
411        remaining_zero_components++;
412        for(short i=current_position;i<_kernel_dimension;i++)
413          if(H[i][j]==BigInt(0))
414            M[i]++;
415      }
416
417    // determine minimal number of such components
418    min=remaining_zero_components;
419    // this is the number of zero components in H[0] and an upper bound
420    // for the M[i]
421    for(short i=current_position;i<_kernel_dimension;i++)
422      if(M[i]<min)
423        min=M[i];
424
425    if(min==(const BigInt&)remaining_zero_components)
426      // all zero components in H[0] are zero in each remaining vector
427      // => desired vector does not exist
428      return 0;
429
430    // add the square of the norm to the vectors with the least common zero
431    // components
432    // discard the others
433    for(short i=current_position;i<_kernel_dimension;i++)
434      if(M[i]!=min)
435        M[i]=-1;
436      else
437        for(short j=0;j<columns;j++)
438          M[i]+=H[i][j]*H[i][j];
439    //  Again, at least one M[i] is positive!
440
441    // determine vector to proceed with
442    // This is the vector with the least common zero components with respect
443    // to H[0], but the smallest possible norm.
444    short min_index=0;
445    for(short i=current_position;i<_kernel_dimension;i++)
446      if(M[i]>BigInt(0))
447        if(min_index==0)
448          min_index=i;
449        else
450          if(M[i]<M[min_index])
451            min_index=i;
452
453    // Now, a multiple of H[min_index] will be added to the already constructed
454    // vector H[0].
455    // For a better handling, it is swapped with the vector at current_position
456    // (only pointers).
457
458    if(min_index!=current_position)
459    {
460      BigInt* swap=H[min_index];
461      H[min_index]=H[current_position];
462      H[current_position]=swap;
463    }
464
465
466// STEP 4: Choose a convenient multiple of H[current_position] to add to H[0].
467// The number of factors "mult" that have to be tested is bounded by the
468// number of nonzero components in H[0] (for each such components, there is at
469// most one such factor that will eliminate it in the linear combination
470// H[0] + mult*H[current_position].
471
472    found=FALSE;
473
474    for(short mult=1;found==FALSE;mult++)
475    {
476      found=TRUE;
477
478      // check if any component !=0 of H[0] becomes zero by adding
479      // mult*H[current_position]
480      for(short j=0;j<columns;j++)
481        if(H[0][j]!=BigInt(0))
482          if(H[0][j]+(const BigInt&)mult*H[current_position][j]
483            ==BigInt(0))
484            found=FALSE;
485
486      if(found==TRUE)
487        for(short j=0;j<columns;j++)
488          H[0][j]+=(const BigInt&)mult*H[current_position][j];
489      else
490        // try -mult
491      {
492
493        found=TRUE;
494
495        // check if any component !=0 of H[0] becomes zero by subtracting
496        // mult*H[current_position]
497        for(short j=0;j<columns;j++)
498          if(H[0][j]!=BigInt(0))
499            if(H[0][j]-(const BigInt&)mult*H[current_position][j]
500              ==BigInt(0))
501              found=FALSE;
502
503        if(found==TRUE)
504          for(short j=0;j<columns;j++)
505            H[0][j]-=(const BigInt&)mult*H[current_position][j];
506      }
507    }
508  }
509
510// When reaching this line, an error must have occurred.
511  cerr<<"FATAL ERROR in short matrix::compute_nonnegative_vector()"<<endl;
512  abort();
513}
514
515short matrix::compute_flip_variables(short*& F)
516{
517  // first compute nonzero vector
518  int okay=compute_nonzero_kernel_vector();
519
520  if(!okay)
521  {
522    cout<<"\nWARNING: short matrix::compute_flip_variables(short*&):\n"
523      "kernel of the matrix contains no vector with nonzero components,\n"
524      "no flip variables computed"<<endl;
525    return -1;
526  }
527
528  // compute variables to flip; these might either be those corresponding
529  // to the positive components of the kernel vector without zero components
530  // or those corresponding to the negative ones
531
532  short r=0;
533  // number of flip variables
534
535  for(short j=0;j<columns;j++)
536    if(H[0][j]<BigInt(0))
537      r++;
538  // remember that all components of H[0] are !=0
539
540  if(r==0)
541    // no flip variables
542    return 0;
543
544  if(2*r>columns)
545    // more negative than positive components in H[0]
546    // all variables corresponding to positive components will be flipped
547  {
548    r=columns-r;
549    F=new short[r];
550    memset(F,0,r*sizeof(short));
551    short counter=0;
552    for(short j=0;j<columns;j++)
553      if(H[0][j]> BigInt(0))
554      {
555        F[counter]=j;
556        counter++;
557      }
558  }
559  else
560    // more (or as many) positive than negative components in v
561    // all variables corresponding to negative components will be flipped
562  {
563    F=new short[r];
564    memset(F,0,r*sizeof(short));
565    short counter=0;
566    for(short j=0;j<columns;j++)
567      if(H[0][j]< BigInt(0))
568      {
569        F[counter]=j;
570        counter++;
571      }
572  }
573
574  return r;
575}
576
577
578
579
580short matrix::hosten_shapiro(short*& sat_var)
581{
582
583  if(_kernel_dimension==-2)
584    // lattice basis not yet computed
585    LLL_kernel_basis();
586
587  if(_kernel_dimension==-1)
588  {
589    cerr<<"\nWARNING: short matrix::hosten_shapiro(short*&):\n"
590      "error in kernel basis, cannot compute the saturation variables"<<endl;
591    return 0;
592  }
593
594  if(_kernel_dimension==0)
595    // the toric ideal corresponding to the kernel lattice is the zero ideal,
596    // no saturation variables necessary
597    return 0;
598
599  // Now, the kernel dimension is positive.
600
601  if(columns==1)
602    // matrix consists of one zero column, kernel is generated by the vector
603    // (1) corresponding to the toric ideal <x-1> which is already staurated
604    return 0;
605
606  short number_of_sat_var=0;
607  sat_var=new short[columns/2];
608  memset(sat_var,0,sizeof(short)*(columns/2));
609
610  BOOLEAN* ideal_saturated_by_var=new BOOLEAN[columns];
611  // auxiliary array used to remember by which variables the ideal has still to
612  // be saturated
613  for(short j=0;j<columns;j++)
614    ideal_saturated_by_var[j]=FALSE;
615
616  for(short k=0;k<_kernel_dimension;k++)
617  {
618    // determine number of positive and negative components in H[k]
619    // corresponding to variables by which the ideal has still to be saturated
620    short pos_sat_var=0;
621    short neg_sat_var=0;
622
623    for(short j=0;j<columns;j++)
624    {
625      if(ideal_saturated_by_var[j]==FALSE)
626      {
627        if(H[k][j]> BigInt(0))
628          pos_sat_var++;
629        else
630          if(H[k][j]< BigInt(0))
631            neg_sat_var++;
632      }
633    }
634
635
636    // now add the smaller set to the saturation variables
637    if(pos_sat_var<=neg_sat_var)
638    {
639      for(short j=0;j<columns;j++)
640        if(ideal_saturated_by_var[j]==FALSE)
641          if(H[k][j]> BigInt(0))
642            // ideal has to be saturated by the variables corresponding
643            // to positive components
644          {
645            sat_var[number_of_sat_var]=j;
646            ideal_saturated_by_var[j]=TRUE;
647            number_of_sat_var++;
648          }
649          else
650            if(H[k][j]< BigInt(0))
651              // then the ideal is automatically saturated by the variables
652              // corresponding to negative components
653              ideal_saturated_by_var[j]=TRUE;
654    }
655    else
656    {
657      for(short j=0;j<columns;j++)
658        if(ideal_saturated_by_var[j]==FALSE)
659          if(H[k][j]< BigInt(0))
660            // ideal has to be saturated by the variables corresponding
661            // to negative components
662          {
663            sat_var[number_of_sat_var]=j;
664            ideal_saturated_by_var[j]=TRUE;
665            number_of_sat_var++;
666          }
667          else
668            if(H[k][j]> BigInt(0))
669              // then the ideal is automatically saturated by the variables
670              // corresponding to positive components
671              ideal_saturated_by_var[j]=TRUE;
672    }
673  }
674
675  // clean up memory
676  delete[] ideal_saturated_by_var;
677
678  return number_of_sat_var;
679}
680
681
682
683
684//////////////////// output ///////////////////////////////////////////////
685
686
687
688
689void matrix::print() const
690{
691  printf("\n%3d x %3d\n",rows,columns);
692
693  for(short i=0;i<rows;i++)
694  {
695    for(short j=0;j<columns;j++)
696      printf("%6d",coefficients[i][j]);
697    printf("\n");
698  }
699}
700
701
702void matrix::print(FILE *output) const
703{
704  fprintf(output,"\n%3d x %3d\n",rows,columns);
705
706  for(short i=0;i<rows;i++)
707  {
708    for(short j=0;j<columns;j++)
709      fprintf(output,"%6d",coefficients[i][j]);
710    fprintf(output,"\n");
711  }
712}
713
714
715void matrix::print(ofstream& output) const
716{
717  output<<endl<<setw(3)<<rows<<" x "<<setw(3)<<columns<<endl;
718
719  for(short i=0;i<rows;i++)
720  {
721    for(short j=0;j<columns;j++)
722      output<<setw(6)<<coefficients[i][j];
723    output<<endl;
724  }
725}
726#endif  // matrix.cc
Note: See TracBrowser for help on using the repository browser.