source: git/IntegerProgramming/matrix.cc @ 6ba162

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