fork download 
  1. #include <iostream>
  2. #include <fstream>
  3. #include <sstream>
  4. #include <vector>
  5. #include <string>
  6. #include <cstdlib>
  7. #include <cassert>
  8. #include <cmath>
  9.  
  10. enum ErrorCode { None = 0, Singular, BadFormat, NotFound };
  11. class Exception {
  12.   enum ErrorCode ErrorCode;
  13. public:
  14.   Exception() : ErrorCode(None) {}
  15.   Exception(enum ErrorCode e) : ErrorCode(e) {}
  16.   enum ErrorCode getErrorCode() { return ErrorCode; }
  17. };
  18.  
  19. class Matrix {
  20. public:
  21.   int n;
  22.   std::vector<std::vector<double> > m;
  23.   Matrix() : n(0) { m.clear(); };
  24.   Matrix(int n, double diag);
  25.   void readFromFile(char const *filename);
  26.  
  27.   friend std::ostream &operator<<(std::ostream &stream, Matrix m);
  28.   friend Matrix operator*(Matrix &a, Matrix &b);
  29. };
  30.  
  31. Matrix::Matrix(int n, double diag) {
  32.   this->n = n;
  33.   std::vector<double> v;
  34.   for (int i = 0; i < n; i++) {
  35.     for (int j = 0; j < n; j++)
  36.       if (i == j)
  37.         v.push_back(diag);
  38.       else
  39.         v.push_back(0.0);
  40.     m.push_back(v);
  41.     v.clear();
  42.   }
  43. }
  44.  
  45. void Matrix::readFromFile(char const *filename) {
  46.   std::ifstream fin(filename);
  47.   if(fin) {
  48.     std::vector<double> line;
  49.     std::string buff;
  50.     int index = 0;
  51.     while (getline(fin, buff)) {
  52.       double data;
  53.       std::stringstream ss;
  54.       ss << buff;
  55.       int i;
  56.       for (i = 0; ss >> data; i++) {
  57.         line.push_back(data);
  58.         ss.ignore();
  59.       }
  60.       m.push_back(line);
  61.       line.clear();
  62.       if (index == 0) {
  63.         this->n = i;
  64.       } else if (i != this->n) {
  65.           std::cerr << "line: " << index << "'s colomn num is not " << this->n << " but " << i << std::endl;
  66.           throw Exception();
  67.       } else if (index >= this->n) {
  68.           std::cerr << "line: " << index << " is too large: " << this->n << std::endl;
  69.           throw Exception();
  70.       }
  71.       index++;
  72.     }
  73.     fin.close();
  74.   } else {
  75.     throw Exception(NotFound);
  76.   }
  77. }
  78.  
  79. std::ostream &operator<<(std::ostream &stream, Matrix m) {
  80.   std::vector<std::vector<double> >::iterator p;
  81.   std::vector<double>::iterator q;
  82.   for (p = m.m.begin(); p != m.m.end(); p++) {
  83.     for (q = (*p).begin(); q != (*p).end(); q++)
  84.       stream << *q << ", ";
  85.     stream << std::endl;
  86.   }
  87.   return stream;
  88. }
  89.  
  90. Matrix operator*(Matrix &a, Matrix &b) {
  91.   Matrix c(a.n, 0.0);
  92.   for (int i = 0; i < a.n; i++)
  93.     for (int j = 0; j < b.n; j++) {
  94.       double s = 0.0;
  95.       for (int k = 0; k < c.n; k++)
  96.         s += a.m[i][k] * b.m[k][j];
  97.       c.m[i][j] = s;
  98.     }
  99.   return c;
  100. }
  101.  
  102. class GaussianEliminationPartial : public Matrix {
  103. public:
  104.   void operator()(Matrix &m, Matrix &e) {
  105.     assert(m.n == e.n);
  106.     /* forward */
  107.     for (int i = 0; i < m.n; i++) {
  108.       /* partial-pivotting */
  109.       double max = std::fabs(m.m[i][i]);
  110.       int k = i;
  111.       for (int j = i + 1; j < m.n; j++) {
  112.         double tmp;
  113.         if (max < (tmp = std::fabs(m.m[j][i]))) {
  114.           max = tmp;
  115.           k = j;
  116.         }
  117.       }
  118.       if (k != i)
  119.         for (int j = 0; j < m.n; j++) {
  120.           double r;
  121.           r = m.m[k][j];
  122.           m.m[k][j] = m.m[i][j];
  123.           m.m[i][j] = r;
  124.           r = e.m[k][j];
  125.           e.m[k][j] = e.m[i][j];
  126.           e.m[i][j] = r;
  127.         }
  128.       /* pivotting finished */
  129.       if (std::fabs(max) < 1.0e-6) {
  130.         throw Exception(Singular);
  131.       }
  132.  
  133.       double r;
  134.       r = m.m[i][i];
  135.       for (int j = 0; j < m.n; j++) {
  136.         m.m[i][j] /= r;
  137.         e.m[i][j] /= r;
  138.       }
  139.       for (int j = i + 1; j < m.n; j++) {
  140.         if ((r = m.m[j][i]) != 0) {
  141.           assert(std::fabs(r) > 0);
  142.           for (int k = 0; k < m.n; k++) {
  143.             m.m[j][k] /= r;
  144.             m.m[j][k] -= m.m[i][k];
  145.             e.m[j][k] /= r;
  146.             e.m[j][k] -= e.m[i][k];
  147.           }
  148.         }
  149.       }
  150.     }
  151.     /* backward */
  152.     for (int i = m.n - 1; i >= 0; --i) {
  153.       for (int j = i - 1; j >= 0; --j) {
  154.         e.m[j][0] -= e.m[i][0] * m.m[j][i];
  155.         for (int k = 1; k < e.n; k++) {
  156.           e.m[j][k] -= e.m[i][k] * m.m[j][i];
  157.         }
  158.         m.m[j][i] = 0.0;
  159.       }
  160.     }
  161.   }
  162. };
  163.  
  164. Matrix inverse(Matrix m) {
  165.   Matrix e(m.n, 1.0);
  166.   GaussianEliminationPartial ge;
  167.   ge(/*ref*/m, /*ref*/e);
  168.   std::cout << m << '\n' << e << std::endl;
  169.   return e;
  170. }
  171.  
  172. char const *fileName = "matrix.csv";
  173. int main() {
  174.   Matrix a, b;
  175.   try {
  176.     a.readFromFile(fileName);
  177.     std::cout << "a = \n" << a << std::endl;
  178.     b = inverse(a);
  179.     std::cout << "b = \n" << b << std::endl;
  180.     std::cout << "a * b = \n" << a * b << std::endl;
  181.   } catch (Exception e) {
  182.     switch (e.getErrorCode()) {
  183.     case NotFound:
  184.       std::cerr << "cannnot open the file." << std::endl;
  185.       break;
  186.     case BadFormat:
  187.       std::cerr << "matrix: bad format." << std::endl;
  188.       break;
  189.     case Singular:
  190.       std::cout << "not have an inverse matrix." << std::endl;    
  191.       break;
  192.     default:
  193.       std::cout << "Internal Error." << std::endl;          
  194.       ;
  195.     }
  196.   }
  197.   return 0;
  198. }
  199. /* end */
  200.  
Success #stdin #stdout #stderr 0s 3040KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
stderr 
cannnot open the file.
https://ideone.com/WSyf7a
language:
C++ (gcc 8.3)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!