CodeArchive LinearAlgebra

From Ubcacm
Jump to: navigation, search

Linear Algebra

Should support Gaussian Elimination, finding determinant, feasibility testing, maximization LP, constructing dual, etc.

gaussian(sys) : Use Gaussian Jordan elimination to solve a system det(M) : Use Gaussian Jordan elimination to find the determinant gauss_sol(M, x) : returns a vector of boolean denoting free variables. In x, return the assignment if all free variables are 0.

create_system(coef, cons, obj) : creates a system for maximization create_dual(coef, cons, obj) : creates a dual system for maximization (the objective is negated) feasible(sys) : converts a LP system into a feasible system (if possible). Return true if successful. get_sol(sys, x) : returns the optimal objective value and the assignment in x (only x.size() elements are returned) solve_from_dual(coef, cons, obj, x) : returns the optimal value and puts the assignment to x. This solves the dual and then recover the primal. This is for educational purpose.

All LP tested with 10498... 

#include <iostream>
#include <vector>
#include <string>
#include <cmath>
using namespace std;

typedef long long ll;
ll gcd(ll a, ll b) {
  if (a < 0) a = -a; if (b < 0) b = -b;
  while(a && b) {
    if (a > b) a%=b;
    else b %= a;
  }
  return a+b;
}

struct rational_t {
  ll nu, de;
  rational_t(const ll &n = 0, const ll &d = 1): nu(n), de(d) {
    ll g = gcd(n, d);
    nu /= g; de /= g;
    if (d < 0) { nu = -nu; de = -de; }
  }
  rational_t operator+(const rational_t& b) const { return rational_t( nu*b.de+de*b.nu, de*b.de ); }
  rational_t operator-(const rational_t& b) const { return rational_t( nu*b.de-de*b.nu, de*b.de ); }
  rational_t operator*(const rational_t& b) const { return rational_t( nu*b.nu, de*b.de ); }
  rational_t operator/(const rational_t& b) const { return rational_t( nu*b.de, de*b.nu ); }
  bool operator == (const rational_t & b) const { return nu * b.de == b.nu * de; }
  bool operator == (const int &k) const { return nu == k * de; }
  bool operator < (const rational_t& b) const { return nu * b.de < b.nu * de; }
  bool operator > (const rational_t& b) const { return nu * b.de > b.nu * de; }
};
rational_t operator-(const rational_t& a) { return rational_t(-a.nu, a.de); }

typedef long double ld;
typedef vector < ld > vd;
typedef vector < vd > vvd;
const ld EPS = 1e-9;
#define ISZERO(x) (fabs(x) < EPS)
#define ISPOS(x) ( x > EPS )
//#define ISZERO(x) ((x).nu == 0)
//#define ISPOS(x) (x.nu > 0)

void pivot(vvd& A, int r, int c) { // A is A|b
  int m = A.size(), n = A[0].size();
  ld tmp = A[r][c];
  for (int i = 0; i < n; ++i) A[r][i] = A[r][i] / tmp;
  for (int i = 0; i < m; ++i)
    if (i != r) {
      ld k = A[i][c];
      for(int j = 0; j < n; ++j)
        A[i][j] = A[i][j] - A[r][j]*k;
    }
}

bool gaussian(vvd& A) {
  int m = A.size(), n = A[0].size() - 1;
  vector<bool> used(m, 0);
  for (int c = 0; c < n; ++c)
    for (int r = 0; r < m; ++r)
      if (!used[r] && !ISZERO(A[r][c]) ) {pivot(A, r, c); used[r] = true; }
  for (int r = 0; r < m; ++r) if(!used[r] && !ISZERO(A[r].back()) ) return false;
  return true;
}

ld det(vvd& A) { // Find the determinant, will destroy matrix...
  if (A.size() != A[0].size()) return 0;
  int n = A.size(), j;
  ld ret = 1;
  for (int i = 0; i < n; ++i) {
    for (j = i; j < n; ++j) if (!ISZERO( A[j][i] )) {
        if (j != i) { A[j].swap(A[i]); ret = ret * -1; }
        ret = ret * A[i][i]; pivot(A, i, i);
        break;
      }
    if (j == n) return 0;
  }
  return ret;
}

vvd create_system(const vvd& coef, const vd &cons, const vd& obj) {
  int m = coef.size(), n = obj.size();
  vvd A(m + 1, vd(n + m + 1, 0));
  for(int i = 0; i < m; i++) for(int j = 0; j < n; j++) A[i][j] = coef[i][j];
  for(int i = 0; i < n; i++) A.back()[i] = -obj[i];
  for(int i = 0; i < m; i++) A[i][n + i] = 1;
  for(int i = 0; i < m; i++) A[i].back() = cons[i];
  return A;
}

vvd create_dual(const vvd& coef, const vd &cons, const vd& obj) {
  int m = obj.size(), n = cons.size();
  vvd A(m + 1, vd(n + m + 1, 0));
  for(int i = 0; i < n; i++) for(int j = 0; j < m; j++) A[j][i] = -coef[i][j];
  for(int i = 0; i < n; i++) A.back()[i] = cons[i]; // Duals have negated objective anyway
  for(int i = 0; i < m; i++) A[i][n + i] = 1;
  for(int i = 0; i < m; i++) A[i].back() = -obj[i];
  return A;
}

bool maximize(vvd& A, int skip = 0){
  int m = A.size(), n = A[0].size();

  while(true){
    int c = -1; ld mn = 0;
    for(int i = 0; i + 1 < n; i++)
      if(A.back()[i] < mn) mn = A.back()[c=i];

    if( ISZERO(mn) ) return true;;

    int r = -1; mn = 1e99; // HUGE
    for(int i = 0; i + 1 + skip < m; i++)
      if( ISPOS(A[i][c]) && A[i].back()/A[i][c] < mn)
        mn = A[i].back()/A[r=i][c];
    if(mn > 1e98) return false; // Still HUGE, Unbounded
    pivot(A, r, c);
  }
  return false;
}

bool feasible(vvd & A) {
  int m = A.size(), n = A[0].size();
  ld tmp;
  for (int i = 0; i < m; ++i) {
    tmp = A[i].back();
    A[i].back() = -1;
    A[i].push_back(tmp);
  }
  A.push_back(vd(n+1, 0));
  A.back()[n-1] = 1; // Objective
  int r = -1; ld mn = 0;
  for (int i = 0; i < m; ++i)
    if (A[i].back() < mn) mn = A[r=i].back();
  if ( ISPOS(-mn) ) {
    pivot(A, r, n-1);
    maximize(A, 1);
    if ( ISPOS(A.back().back()) ) return false;
  }
  A.pop_back();
  for (int i = 0; i < m; ++i) {
    tmp = A[i].back();
    A[i].pop_back();
    A[i].back() = tmp;
  }
  return true;
}

ld get_sol(const vvd& coef, vd& x) { // Put the computed value to x
  for (int i = 0; i < x.size(); ++i) {
    int cnt = 0;
    for (int j = 0; j < coef.size(); ++j)
      if ( !ISZERO(coef[j][i]) ) ++cnt;
    if (cnt != 1) continue;

    for (int j = 0; j < coef.size(); ++j)
      if ( ISZERO(coef[j][i] - 1) ) {
        x[i] = coef[j].back();
        break;
      }
  }
  return coef.back().back();
}

vector<bool> gauss_sol(const vvd& coef, vd&x) { // Return the assignment to each x. Free variables are always 0
  int m = coef.size(), n = coef[0].size()-1;
  vector<bool> ret(n, 0);
  for (int i = 0; i < n; ++i) {
    x[i]=0;
    if (ret[i]) continue;
    int cnt = 0;
    for (int j = 0; j < m; ++j) // Exactly one 1
      if ( !ISZERO(coef[j][i]) ) ++cnt;
    if (cnt != 1) { ret[i] = true; continue; }
    for (int j = 0; j < m; ++j)
      if (ISZERO( coef[j][i] - 1 ) ) {
        x[i] = coef[j].back();
        for (int k = i+1; k < n; ++k)
          if (!ISZERO( coef[j][k] )) ret[k] = true;
        break;
      }
  }
  return ret;
}

ld solve_from_dual(const vvd& coef, const vd& cons, const vd& obj, vd& x) {
  int m = cons.size(), n = obj.size();
  vvd sys(create_dual(coef, cons, obj));
  if(!feasible(sys)) return 1e99; // Unbounded
  maximize(sys); // Solve the dual
  vd y(cons.size(), 0);
  ld opt = -get_sol(sys, y); // We have optimal value and y
  vvd lin;
  for (int i = 0; i < m; ++i) if (!ISZERO(y[i])) {
    lin.push_back(coef[i]); lin.back().push_back(cons[i]); }
  gaussian(lin);
  vector<bool> tmp(gauss_sol(lin, x));
  return opt;
}

-- Main.MatthewChan - 05 Apr 2006


Self-contained determinant: 

ld data[64][64], input[64][64];
char help[64][64];

ld det(int n){
    ld ans = 1;
    for(int j = 0; j < n; j++){
   // find non-zero
   for(int i = j; i < n; i++) if(data[i][j]){
       for(int k = j; k < n; k++) swap(data[i][k], data[j][k]);
       if(i != j) ans *= -1;
       break;
   }
   if(data[j][j]) for(int i = j + 1; i < n; i++){
       ld mult = data[i][j]/data[j][j]; data[i][j] = 0;
       for(int k = j + 1; k < n; k++) data[i][k] -= mult*data[j][k];
   }
    }
    for(int i = 0; i < n; i++) ans *= data[i][i];
    return ans;
}
Retrieved from "http://ubcacm.stuffedcow.net/index.php?title=CodeArchive_LinearAlgebra&oldid=1520"