Re: target/5505: Doubts about a patch for OSF

["Richard B. Kreckel" <kreckel at ginac dot de>]

Hi,

On Fri, 22 Feb 2002, Rainer Orth wrote:
> Richard B. Kreckel writes:
> 
> > Rainer, you probably know by now that this whole issue has nothing to do
> > at all with CLN?  It seems like your patch breaks anything when using
> 
> most likely, yes.
> 
> > -fno-exceptions.  Attached is a trial program that only uses doubles and
> 
> Only in combination with -O2, it seems.
> 
> > some STL containers.  (Sorry for not boiling it down from 100 to 10 lines,
> > I have unfortunately very little time right now.)
> 
> Please try to do so: my knowledge of C++ is close to zero, so a minimal
> testcase would help enormously.

Sure, I will, next week.  Sorry, again I am only for a couple of minutes
in my office...

> > Can you confirm this?  Otherwise we are back with the theory that my boxen
> > are broken / patched wrongly...
> 
> Not yet: I get a SEGV even if compiling without any special options:
> 
> > g++ -o minor minor.cc
> > ./minor
> permanent of 6x6-matrix
> c==4Segmentation fault
> 
> The crash happens in minor.c, l. 80:
> 
> 80                          Pkey[j] = Pkey[j-1]+1;
> 
> j = 0 at this point, so Pkey[j-1] is out of bounds.  The crash happens when
> compiling with Compaq cxx as well.

Sorry, silly me.  I've attached an old case and that index out of bound
has been fixed in the real sources ages ago.  However, the crash I see
seems to have nothing to do with this line!  Modifying the condition such
that j==0 does not happen still results in a crashing program when
compiled with -O2 -fno-exceptions and a working program when compiled with
-O2 alone.  There is an example modified accordingly attached to this
email and this time I have made sure it also runs with -lefence and stuff.

If you have some time, it would be nice if you could confirm the
dependence on -fno-exceptions before I start stripping that down.  Thanks!

Regards
    -richy.
-- 
Richard B. Kreckel
<Richard.Kreckel@Uni-Mainz.DE>
<http://wwwthep.physik.uni-mainz.de/~kreckel/>

#include <cmath>
#include <iostream>
#include <vector>
#include <map>
#include <algorithm>
using namespace std;

int size = 7;

void init(vector<double> & v)
{
    v.clear();
    for (int i=0; i<size*size; ++i)
        v.push_back((10.0*rand())/(RAND_MAX+1.0)-5.0);
}

double det_minor(const vector<double> &v)
{
    // for small matrices the algorithm does not make sense:
    if (size==1)
        return v[0];
    if (size==2)
        return v[0]*v[3]-v[2]*v[1];
    if (size==3)
        return ((v[4]*v[8]-v[5]*v[7])*v[0]-
                (v[1]*v[8]-v[2]*v[7])*v[3]+
                (v[1]*v[5]-v[4]*v[2])*v[6]);

    // we store our subminors in these containers
    typedef map<vector<unsigned>,double> Rmap;
    typedef map<vector<unsigned>,double>::value_type Rmap_value;
    Rmap A, B;
    double det = 0.0;
    vector<unsigned> Pkey;    // Unique flipper counter for the partition
    Pkey.reserve(size);
    vector<unsigned> Mkey;    // key for minor determinant (a partition of Pkey)
    Mkey.reserve(size-1);
    // initialize A with last column:
    for (unsigned r=0; r<size; ++r) {
        Pkey.erase(Pkey.begin(),Pkey.end());
        Pkey.push_back(r);
        A.insert(Rmap_value(Pkey,v[size*r+size-1]));
    }
    // clog << "loop: " << endl;
    for (int c=size-2; c>=0; --c) {
      clog << "c==" << c  << flush;
        Pkey.erase(Pkey.begin(),Pkey.end());  // don't change capacity
        Mkey.erase(Mkey.begin(),Mkey.end());
        for (unsigned i=0; i<size-c; ++i)
            Pkey.push_back(i);
        unsigned fc = 0;  // controls logic for our strange flipper counter
        do {
            det = 0.0;
            for (unsigned r=0; r<size-c; ++r) {
                // maybe there is nothing to do?
                if (v[Pkey[r]*size+c]==0.0)
                    continue;
                // create the sorted key for all possible minors
                Mkey.erase(Mkey.begin(),Mkey.end());
                for (unsigned i=0; i<size-c; ++i)
                    if (i!=r)
                        Mkey.push_back(Pkey[i]);
                // fetch the minors and compute the new determinant
                if (r%2)
                    det -= v[Pkey[r]*size+c]*A[Mkey];
                else
                    det += v[Pkey[r]*size+c]*A[Mkey];
            }
            // Store the new determinant at its place in B:
            if (det!=0.0)
                B.insert(Rmap_value(Pkey,det));
            // increment our strange flipper counter
            for (fc=size-c; fc>0; --fc) {
                ++Pkey[fc-1];
                if (Pkey[fc-1]<fc+c)
                    break;
            }
            if (fc<size-c && fc>0)
                for (unsigned j=fc; j<size-c; ++j)
                    Pkey[j] = Pkey[j-1]+1;
        } while(fc);
        // change the role of A and B:
        A = B;
        B.clear();
	clog << endl;
    }

    return det;
}

int main(void)
{
    srand((unsigned)time(NULL));
    double d0;
    vector<double> m;
    for (size=6; size<11; ++size) {
        cout << "determinant of " << size << "x" << size << "-matrix" << endl;
        for (int i=0; i<20; ++i) {
            init(m);
            d0 = det_minor(m);
            cout << d0 << endl;
        }
    }
}