splitdiag_genpermute.cpp

Go to the documentation of this file.

#include <iostream>
#include <algorithm>
#include <vector>
#include <numeric>
#include <fstream>

using namespace std;
#define PEPS 0.00001
#define NEPS -0.00001

template<typename _ForwardIter, typename T>
static void iota(_ForwardIter __first, _ForwardIter __last, T __val)
{
        while (__first != __last)
                *__first++ = __val++;
}


int main(int argc, char* argv[])
{
        if(argc < 2)
        {
                cout << "Don't forget to pass the input file name" << endl;
        }       
        ifstream input(argv[1]);
        string rhsmatrix = string("restrict_T_") + string(argv[1]);
        string diagonal = string("diag_") + string(argv[1]);
        string nondiagmat = string("offdiag_") + string(argv[1]);

        ofstream rhsout(rhsmatrix.c_str());
        ofstream diagout(diagonal.c_str());
        ofstream matout(nondiagmat.c_str());

        int m, n, nnz;
        input >> m >> n >> nnz;
        vector<int> perm(n);
        iota(perm.begin(), perm.end(), 1);      // 1-based
        random_shuffle(perm.begin(), perm.end());
        int half = m/2;

        if(half *2 < m) // odd
                rhsout << m << "\t" << (half+1) << "\t" << m << endl;
        else    // even
                rhsout << m << "\t" << half << "\t" << m << endl;
        
        int one = 1;
        for(int i=0; i< half; ++i)
        {
                rhsout << perm[2*i] << "\t" << (i+1) << "\t" << one << endl;
                rhsout << perm[2*i+1] << "\t" << (i+1) << "\t" << one << endl;  
        }       
        if (half *2 < m)        // m was odd
        {
                rhsout << perm.back() << "\t" << (half+1) << "\t" << one << endl;
        }
        int row, column;
        double value;
        int diagnnz = 0;
        int nondiagnnz = 0;
        int zero = 0;
        diagout << m << "\t" << one << "\t" << zero << endl;
        matout << m << "\t" << n << "\t" << zero << endl;
        while(!input.eof())
        {
                input >> row >> column >> value;
                // keep explicit zeros
                if(row == column)
                {
                        diagout << row << "\t" << one << "\t" << value << endl;
                        diagnnz++;
                }
                else
                {
                        matout << row << "\t" << column << "\t" << value << endl;
                        nondiagnnz++;
                }
        }
        cout << "Diagonal nnz: " << diagnnz << endl;
        cout << "Off-diagonal nnz: " << nondiagnnz << endl;
        // will need to reopen this and write first lines
        rhsout.close();
        diagout.close();
        matout.close();

        return 0;
}