fsi_output.cpp

#include "pugixml.hpp"
#include <armadillo>  // linear algebra library
#include <cassert>

#if 0
#include <vtkNew.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkStructuredGrid.h>
#include <vtkType.h>
#include <vtkUnstructuredGrid.h>
#include <vtkStructuredGridWriter.h>
#include <vtkUnstructuredGridWriter.h>
#endif

#include "CResults.h"
#include "importExport.h"
#include "DSimpleMessage.h"
#include "fsi_output.hpp"

#ifndef VLM_WRITE_PATH
#ifdef Q_OS_WIN32
#define VLM_WRITE_PATH "./" //TODO remove before release
#else
#define VLM_WRITE_PATH "/home/wm/src/WLM_PM_CODE/data/"
//#define VLM_WRITE_PATH "/home/wm/ml/WLM_PM_CODE/data/"
//#define VLM_WRITE_PATH "/tmp/VLMdata"
#endif
#endif

#ifdef FSI_DEBUG
#define PRECISION "%1.7g"  // printf precision to standardize debugging output
#endif

#define NROWS(matrix) (matrix).n_rows
#define NCOLS(matrix) (matrix).n_cols
#define N3D(matrix) (matrix).n_slices
#define ENDROW(var) (var)((var).n_rows - 1  // note: unusual looks - missing close paren
#define ENDCOL(var) (var).n_cols - 1
#define ENDN3D(var) (var).n_slices - 1

using namespace std;
using namespace arma;

#define SP(var, cls) vtkSmartPointer<cls> var = vtkSmartPointer<cls>::New()
#define COMP(i, j) slice(j).col(i)

#ifndef _WIN32
#include <cerrno>
#include <cctype>
#define fprintf_s fprintf
error_t fopen_s(FILE **f, const char *name, const char *mode) {
    error_t ret = 0;
    assert(f);
    *f = fopen(name, mode);
    /* Can't be sure about 1-to-1 mapping of errno and MS' errno_t */
    if (!*f)
        ret = errno;
    return ret;
}
#endif

/*
 * generate_nastran is for internal use only. It generates a NASTRAN .bdf mesh according to
 * boundary points generated by the solver, which functions as a HAWT solver when run without
 * a NASTRAN file as input.
 *
 * ydist is assumed to already be factored in
 * NOTE: if rotations are selected and factored into the input, they will be left here,
 * although they were not in the original mesh. To use this generator, user must zero rotations
 */
void generate_nastran(const cube &ydistXYZ, const bool include_wake_panel,
                      const double ydist,
                      const double roll,
                      const double AOA,
                      const double sweep,
                      const double span,
                      const double c,
                      vector <string> &comment) {
    // for a 2m + 1 x n + 1 grid, 2m + 2, n + 2
    ///cube XYZ(3, 2 * m + 2, n + 1);  // boundary points
    ///
    cube XYZ(ydistXYZ);
    XYZ.row(1) -= ydist;  // subtract ydist that is assumed to be added in already

    CMesh mesh;
    //SP(og, vtkUnstructuredGrid);
    //SP(wrt, vtkUnstructuredGridWriter);
    //SP(og, vtkStructuredGrid);
    //SP(wrt, vtkStructuredGridWriter);
    //SP(points, vtkPoints);
    //SP(cells, vtkCellArray);
    //Col<vtkIdType> r(4);
    //Cube<vtkIdType> cel(4, XYZ.n_cols - 1, XYZ.n_slices - 1);
    //cel.zeros();
    uword i, j,
          m2 = XYZ.n_cols - 1,  // upper + lower
          n = XYZ.n_slices - 1;
    if (!include_wake_panel) m2--;
    vec::fixed<3> xyz;
    // generate nodes: NASTRAN GRID
    for (j = 0; j < n + 1; j++) {  // last row of points is redundant, since they wrap
        for (i = 0; i < m2; i++) {
            xyz = XYZ.COMP(i, j);
            mesh.insertNodeWithCoords(xyz.memptr());  // GRID
        }
    }
    vecUInt ids(4);
    cerr << "generating nastran, m*2=" << m2 << " n=" << n << endl;
    int ncells = 0;
    // generate elements: NASTRAN CQUAD4
    for (j = 0; j < n; j++) {
        for (i = 0; i < m2; i++) {
            ids[0] = m2 * j + i + 1;  // one-based
            ids[1] = (i < m2 - 1 ? ids[0] + 1 : m2 * j + 1);
            ids[3] = ids[0] + m2;
            ids[2] = (i < m2 - 1 ? ids[3] + 1 : ids[1] + m2);
            //cerr << i << "," << j << ": " << ids[0] << " " << ids[1]<< " " << ids[2]<< " " << ids[3] << endl;
            mesh.insertElemWithNodeIDs(1, ids);  // CQUAD4
            ncells++;
        }
    }
    const string nfn(
#ifdef NASTRAN_GEN_PATH
                NASTRAN_GEN_PATH
#else
#ifdef Q_OS_WIN32
                "c:/Users/Wm/cshared/"
#else
                "/data/samples/"
                //"/media/sf_cshared/"
#endif
#endif
                "sample_mesh_m" + to_string(m2 / 2) + "n" + to_string(n) + ".bdf");
    cerr << "exporting NASTRAN format file: " << nfn << endl;
    comment.emplace(comment.begin(), "# span-wise panels (n) = " + to_string(n));
    comment.emplace(comment.begin(), "# chord-wise panels (m) = " + to_string(m2 / 2));
    int ret = exportNastranFile(nfn,
                                comment,
                                &mesh,
                                nullptr,
                                true,
                                false,
                                false,
                                string(),
                                string(),
                                string(),
                                string(),
                                string());
    cerr << "export nastran returned " << ret << endl;
    cerr << "ncells " << i << " x " << j << " = " << ncells << endl;
}

void structured_grid_out(pugi::xml_node &doc,
                         const char *name,
                         const uword iblade,
                         const cube &bx,
                         const int t,
                         const int parmnrows,
                         bool t_is_name_only) {

    int nrows;

    string buf(name);
    buf += to_string(iblade + 1);  // one-based
    if (t != -1) {
        buf += "_T" + to_string(t + 1);
    }
    if (t == -1 || t_is_name_only) {  // no time step - output entire matrix
        nrows = (parmnrows == -1 ? NROWS(bx) : parmnrows);
        assert(nrows > 0);
    } else {  // end output with indicated time step
        nrows = t + 1;
        if (nrows < 1) {
            cerr << "ERROR: no output for " << buf << endl;
        }
        assert(nrows >= 1);  // must have at least one row to write
        if (nrows > NCOLS(bx)) {
            cerr << "ERROR: not enough input rows (" << NCOLS(bx) <<
                ") to output " << buf << endl;
        }
        assert(nrows <= NCOLS(bx));
    }
    assert(nrows > 0);
    size_t n = 0;
    // TODO: seems that m and n are reversed in extents wrt "3"
    const string extents(to_string(nrows) + " " + to_string(N3D(bx)) + " 3");  // data is 3-tuples
    pugi::xml_node node = doc.append_child(buf.c_str());
    assert(node != nullptr);
    pugi::xml_attribute attr = node.append_attribute("extents");
    attr.set_value(extents.c_str());
    string datastr("\n" INDENT2);  // start at left margin for readability
    for (size_t j = 0; j < N3D(bx); j++) {
        for (size_t i = 0; i < nrows; i++) {
            datastr += to_string(bx(0, i, j)) + " " +
                       to_string(bx(1, i, j)) + " " +
                       to_string(bx(2, i, j)) + " ";
            n += 3;
        }
        datastr += (j + 1 == N3D(bx) ? "\n" INDENT : "\n" INDENT2);
    }
    attr = node.append_attribute("size");
    attr.set_value(n);
    node.append_child(pugi::node_pcdata).set_value(datastr.c_str());
    //cerr << "Wrote grid " << buf << " to xml DOM." << endl;
}

void structured_grid_out(pugi::xml_node &doc,
                         const char *name,
                         const mat &bx, const mat &by, const mat &bz,
                         const int t,
                         const int parmnrows,
                         bool t_is_name_only) {

    int nrows;

    string buf(name);
    if (t != -1) {
        buf += "_T" + to_string(t + 1);
    }
    assert(NROWS(bx) == NROWS(by) && NCOLS(bx) == NCOLS(by));
    assert(NROWS(bx) == NROWS(bz) && NCOLS(bx) == NCOLS(bz));

    if (t == -1 || t_is_name_only) {  // no time step - output entire matrix
        nrows = (parmnrows == -1 ? NROWS(bx) : parmnrows);
        assert(nrows > 0);
    } else {  // end output with indicated time step
        nrows = t + 1;
        if (nrows < 1) {
            cerr << "ERROR: no output for " << buf << endl;
        }
        assert(nrows >= 1);  // must have at least one row to write
        if (nrows > NROWS(bx)) {
            cerr << "ERROR: not enough input rows (" << NROWS(bx) <<
                ") to output " << buf << endl;
        }
        assert(nrows <= NROWS(bx));
    }
    assert(nrows > 0);
    size_t n = 0;
    const string extents(to_string(nrows) + " " + to_string(NCOLS(bx)) + " 3");  // data is 3-tuples
    pugi::xml_node node = doc.append_child(buf.c_str());
    assert(node != nullptr);
    pugi::xml_attribute attr = node.append_attribute("extents");
    attr.set_value(extents.c_str());
    string datastr("\n" INDENT2);  // start at left margin for readability
    for (size_t j = 0; j < NCOLS(bx); j++) {
        for (size_t i = 0; i < nrows; i++) {
            datastr += to_string(bx(i, j)) + " " +
                       to_string(by(i, j)) + " " +
                       to_string(bz(i, j)) + " ";
            n += 3;
        }
        datastr += (j + 1 == NCOLS(bx) ? "\n" INDENT : "\n" INDENT2);
    }
    attr = node.append_attribute("size");
    attr.set_value(n);
    node.append_child(pugi::node_pcdata).set_value(datastr.c_str());
    //cerr << "Wrote grid " << buf << " to xml DOM." << endl;
}

/*!
 * array_data_new(filename, mat x, y, z, optional time step +1)
 *      output to xml DOM arrays X, Y, Z with optional time step
 *      time step t is 1-based, as are output file names
 *      appends to XML doc
 */
void array_out(pugi::xml_node &doc,
               const char *name,
               const vec &bx,
               const int t) {

    int nrows;

    if (t == -1) {  // no time step - output entire matrix
        nrows = NROWS(bx);
        assert(nrows > 0);
    } else {  // end output with indicated time step
        nrows = t + 1;
        if (nrows < 1) {
            cerr << "ERROR: no output for array " << name << endl;
        }
        assert(nrows >= 1);  // must have at least one row to write
        if (nrows > NROWS(bx)) {
            cerr << "ERROR: not enough input rows (" << NROWS(bx) <<
                ") to array " << name << endl;
        }
        assert(nrows <= NROWS(bx));
    }
    string datastr("\n" INDENT2);  // start at left margin for readability
    string extents(to_string(NROWS(bx)));
    pugi::xml_node node = doc.append_child(name);
    assert(node != nullptr);
    pugi::xml_attribute attr = node.append_attribute("extents");
    attr.set_value(extents.c_str());
    for (size_t i = 0; i < nrows; i++) {
        datastr += to_string(bx[i]) + " ";
    }
    datastr += "\n" INDENT;
    attr = node.append_attribute("size");
    attr.set_value(nrows);
    node.append_child(pugi::node_pcdata).set_value(datastr.c_str());
    //cerr << "Wrote master array " << name << " to xml DOM." << endl;
}

/*!
 * matrix_out(xml DOM, varname, mat x, time step)
 *      output to xml DOM file for arrays X, Y, Z with optional time step
 *      assumes that time step is 3rd dimension or matrix
 *      time step t is 1-based TODO ???, as are output file names
 */
void matrix_out(pugi::xml_node &doc,
                const char *name,
                const mat &bx,
                const int t) {

    string buf(name);
    if (t != -1) {
        buf += "_T" + to_string(t + 1);
    }
    const string extents(to_string(NROWS(bx)) + " " + to_string(NCOLS(bx)));
    pugi::xml_node node = doc.append_child(buf.c_str());
    assert(node != nullptr);
    pugi::xml_attribute attr = node.append_attribute("extents");
    attr.set_value(extents.c_str());
    string datastr("\n" INDENT2);  // start at left for readability
    size_t n = 0;
    for (size_t j = 0; j < NCOLS(bx); j++) {
        for (size_t i = 0; i < NROWS(bx); i++) {
            datastr += to_string(bx(i, j)) + " ";
            n++;
        }
        datastr += (j + 1 == NCOLS(bx) ? "\n" INDENT : "\n" INDENT2);
    }
    attr = node.append_attribute("size");
    attr.set_value(n);
    node.append_child(pugi::node_pcdata).set_value(datastr.c_str());
    //cerr << "Wrote matrix " << buf << " to master xml DOM." << endl;
}

/*!
 * matrix_out(filename, cube x, time step)
 *      output to xml DOM arrays X, Y, Z with optional time step
 *      assumes that time step is 3rd dimension or matrix
 *      time step t is 1-based TODO ???, as are output file names
 */
void matrix_out(pugi::xml_node &doc,
                const char *name,
                const cube &bx,
                const int t) {

    size_t n = 0;
    const string extents(to_string(NROWS(bx)) + " " + to_string(NCOLS(bx)));
    pugi::xml_node node = doc.append_child(name);
    assert(node != nullptr);
    pugi::xml_attribute attr = node.append_attribute("extents");
    attr.set_value(extents.c_str());
    string datastr("\n" INDENT2);  // start at left margin for readability
    for (size_t j = 0; j < NCOLS(bx); j++) {
        for (size_t i = 0; i < NROWS(bx); i++) {
            datastr += to_string(bx(i, j, t)) + " ";
            n++;
        }
        datastr += (j + 1 == NCOLS(bx) ? "\n" INDENT : "\n" INDENT2);
    }
    attr = node.append_attribute("size");
    attr.set_value(n);
    node.append_child(pugi::node_pcdata).set_value(datastr.c_str());
    //cerr << "Wrote array " << name << " to xml DOM." << endl;
}

/*!
 * matrix_out(doc, name, matrix_t x, time step)
 *      output to xml DOM arrays
 *      assumes that time step is 1st dimension of matrix
 */
void matrix_out_t(pugi::xml_node &doc,
                  const char *XXfile,
                  const mat &bx,
                  const int t) {

    size_t n = 0;
    string name(XXfile);
    name += "_T" + to_string(t + 1);

    const string extents(to_string(t + 1) + " " + to_string(NCOLS(bx)));
    pugi::xml_node node = doc.append_child(name.c_str());
    assert(node != nullptr);
    pugi::xml_attribute attr = node.append_attribute("extents");
    attr.set_value(extents.c_str());
    string datastr("\n" INDENT2);  // start at left margin for readability
    for (size_t j = 0; j < NCOLS(bx); j++) {
        for (size_t i = 0; i <= t; i++) {
            datastr += to_string(bx(i, j)) + " ";
            n++;
        }
        datastr += (j + 1 == NCOLS(bx) ? "\n" INDENT : "\n" INDENT2);
    }
    attr = node.append_attribute("size");
    attr.set_value(n);
    node.append_child(pugi::node_pcdata).set_value(datastr.c_str());
    //cerr << "Wrote array " << name << " to xml DOM." << endl;
}

#ifdef FSI_DEBUG

void genpath(string &buf, const char *path, const char *XXfile, const int t) {
    buf = path;
    buf += "/";
    buf += XXfile;
    if (t != -1) {
        buf += "_T" + to_string(t + 1);
    }
    buf += ".txt";
}

void genpath_at_t(string &buf, const char *path, const char *XXfile, const int t) {
    buf = path;
    buf += "/";
    buf += XXfile;
    buf += "_at_T" + to_string(t + 1) + ".txt";
}

// TODO this doesn't work right
void outdb(const char *XXfile, const uword xyz, const double *bx, const int m, const int n, const int t) {
    string buf;
    genpath(buf, VLM_WRITE_PATH, XXfile, t);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (size_t i = 0; i < m; i++) {
        for (size_t j = 0; j < n; j++, bx += 3) {
            fprintf(fp, PRECISION, bx[xyz]);
            if (j != n)
                fprintf(fp, " ");
        }
        fprintf(fp, "\n");
    }
    fclose(fp);
    cerr << "Wrote buffer debugging file:" << buf << endl;
}

void outdb(const char *XXfile, const uword xyz, const cube &bx, const int t) {
    string buf;
    genpath(buf, VLM_WRITE_PATH, XXfile, t);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (size_t i = 0; i < NCOLS(bx); i++) {
        for (size_t j = 0; j < N3D(bx); j++) {
            fprintf(fp, PRECISION, bx(xyz, i, j));
            if (j != ENDN3D(bx))
                fprintf(fp, " ");
        }
        fprintf(fp, "\n");
    }
    fclose(fp);
    //cerr << "Wrote cube debugging file:" << buf << endl;
}

void outdb(const char *XXfile, const uword xyz, const mat &bx, const int t) {
    string buf;
    genpath(buf, VLM_WRITE_PATH, XXfile, t);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (size_t j = 0; j < NCOLS(bx); j++) {
        fprintf(fp, PRECISION, bx(xyz, j));
        if (j != ENDCOL(bx))
            fprintf(fp, " ");
    }
    fprintf(fp, "\n");
    fclose(fp);
    //cerr << "Wrote cube debugging file:" << buf << endl;
}

void outdb(const char *XXfile, const uword xyz, const vec &bx, const int t) {
    string buf;
    FILE *fp;
    genpath(buf, VLM_WRITE_PATH, XXfile, t);
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (size_t j = 0; j < NCOLS(bx); j++) {
        fprintf_s(fp, PRECISION, bx(xyz, j));
        if (j != ENDCOL(bx))
            fprintf(fp, " ");
    }
    fprintf_s(fp, "\n");
    fclose(fp);
    cerr << "Wrote vec debugging file:" << buf << endl;
}

void outd(const char *XXfile, const mat &bx, const int t) {
    string buf;
    genpath(buf, VLM_WRITE_PATH, XXfile, t);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (size_t i = 0; i < NROWS(bx); i++) {
        for (size_t j = 0; j < NCOLS(bx); j++) {
            fprintf(fp, PRECISION, bx(i, j));
            if (j != ENDCOL(bx))
                fprintf(fp, " ");
        }
        fprintf(fp, "\n");
    }
    fclose(fp);
    //cerr << "Wrote matrix debugging file:" << buf << endl;
}

void outd(const char *XXfile, const cube &bx, const int t) {
    string buf;
    genpath(buf, VLM_WRITE_PATH, XXfile, t);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (size_t i = 0; i < NROWS(bx); i++) {
        for (size_t j = 0; j < NCOLS(bx); j++) {
            for (size_t k = 0; k < N3D(bx); k++) {
                fprintf(fp, PRECISION, bx(i, j, k));
                if (j != ENDCOL(bx) || k != N3D(bx) - 1)
                    fprintf(fp, " ");
            }
        }
        fprintf(fp, "\n");
    }
    fclose(fp);
    cerr << "Wrote cube debugging file:" << buf << endl;
}

/*!
 * outdatt(file, mat, time t)
 * print 3-D boost matrix at time - time is not a dimension in the matrix
 * for debugging only
 */
void outdatt(const char *XXfile, const cube &bx, const int t) {

    string buf;
    genpath_at_t(buf, VLM_WRITE_PATH, XXfile, t);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (size_t i = 0; i < NROWS(bx); i++) {
        for (size_t j = 0; j < NCOLS(bx); j++) {
            fprintf(fp, PRECISION, bx(i, j, t));
            if (j != ENDCOL(bx))
                fprintf(fp, " ");
        }
        fprintf(fp, "\n");
    }
    fclose(fp);
    //cerr << "Wrote debugging file:" << buf << endl;
}

/*!
 * outdatt(file, mat, time t)
 * print Octave matrix at time with "_AT_" in file name
 * time is not a dimension in the matrix
 * for debugging only - compatibility with MATLAB debugging function
 */
void outdatt(const char *XXfile, const mat &bx, const int t) {

    string buf;
    genpath_at_t(buf, VLM_WRITE_PATH, XXfile, t);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (uword i = 0; i < NROWS(bx); i++) {
        for (uword j = 0; j < NCOLS(bx); j++) {
            fprintf(fp, PRECISION, bx(i, j));
            if (j != NCOLS(bx) - 1)
                fprintf(fp, " ");
        }
        fprintf(fp, "\n");
    }
    fclose(fp);
    //cerr << "Wrote debugging file:" << buf << endl;
}

/*!
 * outd(file, colvec, time optional)
 * print Octave colvec, optionally at time
 * time is not a dimension in the matrix
 * if t provided, append "_T" + time step plus one to file name before extension
 * for debugging only
 */
void outd(const char *XXfile, const vec &bx, const int t) {

    string buf;
    genpath(buf, VLM_WRITE_PATH, XXfile, t);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (uword j = 0; j < NROWS(bx); j++) {
        fprintf(fp, PRECISION, bx(j));
        if (j != NROWS(bx) - 1)
            fprintf(fp, " ");
    }
    fclose(fp);
    //cerr << "Wrote vector debugging file:" << buf << endl;
}

/*!
 * outdcolvec(file, colvec, time optional)
 * print arma colvec, optionally at time
 * - specifically to match output of MATLAB column vectors
 * - no difference in arma vec and colvec
 * time is not a dimension in the matrix
 * if t provided, append "_T" + time step plus one to file name before extension
 * for debugging only
 */
void outdcolvec(const char *XXfile, const colvec &bx, const int t) {

    string buf;
    genpath(buf, VLM_WRITE_PATH, XXfile, t);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (uword j = 0; j < NROWS(bx); j++) {
        fprintf(fp, PRECISION "\n", bx(j));
    }
    fclose(fp);
    //cerr << "Wrote colvec debugging file:" << buf << endl;
}

void outdrowvec(const char *XXfile, const rowvec &bx, const int t) {

    string buf;
    genpath(buf, VLM_WRITE_PATH, XXfile, t);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (uword j = 0; j < bx.size(); j++) {
        fprintf(fp, PRECISION "\n", bx(j));
    }
    fclose(fp);
    //cerr << "Wrote colvec debugging file:" << buf << endl;
}

/*!
 * out2d(file, double *, rows, columns)
 * print 2D array given double *
 * for debugging only
 */
void out2d(const char *XXfile, const double *pbx, const int nrows, const int ncols) {
    string buf;
    genpath(buf, VLM_WRITE_PATH, XXfile);
    FILE *fp;
    if (fopen_s(&fp, buf.c_str(), "w") != 0) {
        cerr << "Cannot write to " << buf << endl;
        return;
    }
    for (int i = 0; i < nrows; i++) {
        for (int j = 0; j < ncols; j++) {
            fprintf(fp, PRECISION, pbx[i * ncols + j]);
            if (j != ncols - 1)
                fprintf(fp, " ");
        }
        fprintf(fp, "\n");
    }
    fclose(fp);
    //cerr << "Wrote debugging file:" << buf << endl;
}

void dumpar(const char *label, const colvec &a) {
    cerr << label << " size=" << a.size() << ":";
    for (unsigned int i = 0; i < a.size(); i++) {
        cerr << a(i) << " ";
    }
    cerr << endl;
}

void dumpmat(const char *label, const mat &a) {
    cerr << label << " size=[" << NROWS(a) << ":" << NCOLS(a) << "]:" << endl;
    for (size_t i = 0; i < NROWS(a); i++) {
        for (size_t j = 0; j < NCOLS(a); j++) {
            printf("%.4f ", a(i, j));
            //cerr << a[i][j] << " ";
        }
        cerr << endl;
    }
}
#endif