source: git/IntegerProgramming/matrix.cc @ c265b9

spielwiese
Last change on this file since c265b9 was 056f98, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes:new +memset git-svn-id: file:///usr/local/Singular/svn/trunk@5377 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 16.9 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    memset(F,0,r*sizeof(short));
555    short counter=0;
556    for(short j=0;j<columns;j++)
557      if(H[0][j]>0)
558      {
559        F[counter]=j;
560        counter++;
561      }
562  }
563  else
564    // more (or as many) positive than negative components in v
565    // all variables corresponding to negative components will be flipped
566  {
567    F=new short[r];
568    memset(F,0,r*sizeof(short));
569    short counter=0;
570    for(short j=0;j<columns;j++)
571      if(H[0][j]<0)
572      {
573        F[counter]=j;
574        counter++;
575      }
576  }
577
578  return r;
579}
580
581
582
583
584short matrix::hosten_shapiro(short*& sat_var)
585{
586
587  if(_kernel_dimension==-2)
588    // lattice basis not yet computed
589    LLL_kernel_basis();
590
591  if(_kernel_dimension==-1)
592  {
593    cerr<<"\nWARNING: short matrix::hosten_shapiro(short*&):\n"
594      "error in kernel basis, cannot compute the saturation variables"<<endl;
595    return 0;
596  }
597
598  if(_kernel_dimension==0)
599    // the toric ideal corresponding to the kernel lattice is the zero ideal,
600    // no saturation variables necessary
601    return 0;
602
603  // Now, the kernel dimension is positive.
604
605  if(columns==1)
606    // matrix consists of one zero column, kernel is generated by the vector
607    // (1) corresponding to the toric ideal <x-1> which is already staurated
608    return 0;
609
610  short number_of_sat_var=0;
611  sat_var=new short[columns/2];
612  memset(sat_var,0,sizeof(short)*(columns/2));
613
614  BOOLEAN* ideal_saturated_by_var=new BOOLEAN[columns];
615  // auxiliary array used to remember by which variables the ideal has still to
616  // be saturated
617  for(short j=0;j<columns;j++)
618    ideal_saturated_by_var[j]=FALSE;
619
620  for(short k=0;k<_kernel_dimension;k++)
621  {
622    // determine number of positive and negative components in H[k]
623    // corresponding to variables by which the ideal has still to be saturated
624    short pos_sat_var=0;
625    short neg_sat_var=0;
626
627    for(short j=0;j<columns;j++)
628    {
629      if(ideal_saturated_by_var[j]==FALSE)
630      {
631        if(H[k][j]>0)
632          pos_sat_var++;
633        else
634          if(H[k][j]<0)
635            neg_sat_var++;
636      }
637    }
638
639
640    // now add the smaller set to the saturation variables
641    if(pos_sat_var<=neg_sat_var)
642    {
643      for(short j=0;j<columns;j++)
644        if(ideal_saturated_by_var[j]==FALSE)
645          if(H[k][j]>0)
646            // ideal has to be saturated by the variables corresponding
647            // to positive components
648          {
649            sat_var[number_of_sat_var]=j;
650            ideal_saturated_by_var[j]=TRUE;
651            number_of_sat_var++;
652          }
653          else
654            if(H[k][j]<0)
655              // then the ideal is automatically saturated by the variables
656              // corresponding to negative components
657              ideal_saturated_by_var[j]=TRUE;
658    }
659    else
660    {
661      for(short j=0;j<columns;j++)
662        if(ideal_saturated_by_var[j]==FALSE)
663          if(H[k][j]<0)
664            // ideal has to be saturated by the variables corresponding
665            // to negative components
666          {
667            sat_var[number_of_sat_var]=j;
668            ideal_saturated_by_var[j]=TRUE;
669            number_of_sat_var++;
670          }
671          else
672            if(H[k][j]>0)
673              // then the ideal is automatically saturated by the variables
674              // corresponding to positive components
675              ideal_saturated_by_var[j]=TRUE;
676    }
677  }
678
679  // clean up memory
680  delete[] ideal_saturated_by_var;
681
682  return number_of_sat_var;
683}
684
685
686
687
688//////////////////// output ///////////////////////////////////////////////
689
690
691
692
693void matrix::print() const
694{
695  printf("\n%3d x %3d\n",rows,columns);
696
697  for(short i=0;i<rows;i++)
698  {
699    for(short j=0;j<columns;j++)
700      printf("%6d",coefficients[i][j]);
701    printf("\n");
702  }
703}
704
705
706void matrix::print(FILE *output) const
707{
708  fprintf(output,"\n%3d x %3d\n",rows,columns);
709
710  for(short i=0;i<rows;i++)
711  {
712    for(short j=0;j<columns;j++)
713      fprintf(output,"%6d",coefficients[i][j]);
714    fprintf(output,"\n");
715  }
716}
717
718
719void matrix::print(ofstream& output) const
720{
721  output<<endl<<setw(3)<<rows<<" x "<<setw(3)<<columns<<endl;
722
723  for(short i=0;i<rows;i++)
724  {
725    for(short j=0;j<columns;j++)
726      output<<setw(6)<<coefficients[i][j];
727    output<<endl;
728  }
729}
730#endif  // matrix.cc
Note: See TracBrowser for help on using the repository browser.