//
// ginac_linalg_funcs.cpp
//
// Warren Weckesser
// Department of Mathematics
// Colgate University
// Hamilton, NY 13346
//

#include <iostream>
#include <ginac/ginac.h>

using namespace std;
using namespace GiNaC;

matrix mminor(matrix A, unsigned r, unsigned c)
    {
    unsigned n = A.rows();
    if (n != A.cols())
        {
        cout << "mminor: matrix is not square.\n";
        exit(-1);
        }
    if (r >= n || c >= n)
        {
        cout << "mminor: row or column out of range.\n";
        exit(-1);
        }
    matrix M = matrix(n-1,n-1);
    for (int i = 0; i < n-1; ++i)
        {
        int i1 = i;
        if (i1 >= r)
            i1 = i+1;
        for (int j = 0; j < n-1; ++j)
            {
            int j1 = j;
            if (j1 >= c)
                j1 = j+1;
            M(i,j) = A(i1,j1);
            }
        }
    return(M);
    }

matrix adjoint(matrix A)
    {
    unsigned n = A.rows();
    if (n != A.cols())
        {
        cout << "adjoint: matrix is not square.\n";
        exit(-1);
        }
    matrix M = matrix(n,n);
    int sr = 1;
    for (int i = 0; i < n; ++i)
        {
        int sc = 1;
        for (int j = 0; j < n; ++j)
            {
            matrix SubM = mminor(A,j,i);
            M(i,j) = sr*sc*SubM.determinant();
            sc = -sc;
            }
         sr = -sr;
         }
    return(M);
    }


//
// jacobian assumes that f is a column vector (n x 1).
//

matrix jacobian(matrix f, lst vars)
    {
    int nrows = f.rows();
    int ncols = vars.nops();
    matrix J = matrix(nrows,ncols);
    for (int i = 0; i < nrows; ++i)
        for (int j = 0; j < ncols; ++j)
            {
            symbol v = ex_to<symbol>(vars[j]);
            J(i,j) = f(i,0).diff(v);
            }
    return J;
    }