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

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

Hi,

On Wed, 20 Feb 2002, I wrote:
[...]
> Err, while trying to debug into the problem I discovered something
> that had escaped my attantion until now: compiling CLN and an example
> (examples/e, or tests/tests or whatever, never mind) with either -O1, -O2
> or -O1 -g resulted in a working test (the linker warnings are of course
> still present) while -O2 -fno-exceptions produced a crashing program.
> Also, these funny warnings:
>  as1: Warning: /tmp/ccb8ZbYD.s, line 6: macro instruction used $at
> appear only when I disable exceptions.  I hadn't noticed it so far because
> I *always* export CXXFLAGS="-O2 -fno-exceptions" prior to building
> CLN.  May I ask you how you configured and tested CLN?  You did not
> specify -fno-exceptions, did you?  Does it work when you do so?

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
-fno-exceptions.  Attached is a trial program that only uses doubles and
some STL containers.  (Sorry for not boiling it down from 100 to 10 lines,
I have unfortunately very little time right now.)

Can you confirm this?  Otherwise we are back with the theory that my boxen
are broken / patched wrongly...

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)
                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 << "permanent of " << size << "x" << size << "-matrix" << endl;
        for (int i=0; i<20; ++i) {
            init(m);
            d0 = det_minor(m);
            cout << d0 << endl;
        }
    }
}