Matrix Operations Calculator (cpp, written by Codex)

Codex wrote this cpp project from the task below. It installed and ran on a clean Ubuntu 22.04 machine as written.

Task given to the agent:

TASK: Matrix Operations Calculator

Write a program that performs common matrix operations including addition, multiplication, transposition, determinant calculation, inversion, and eigenvalue decomposition on matrices loaded from files.

FUNCTIONAL REQUIREMENTS:
- Accept a matrix data file (CSV or JSON format) as a command-line argument
- Support operations selectable via --operation flag: add, subtract, multiply, transpose, determinant, inverse, eigenvalues, rank, trace, and LU decomposition
- For binary operations (add, subtract, multiply), accept a second matrix file via --matrix2 flag
- Support scalar operations: scalar multiplication via --scalar flag applied to the matrix
- Compute matrix properties: dimensions, rank, trace, is-symmetric, is-positive-definite, condition number
- Handle matrices of arbitrary size (up to practical memory limits)
- Support sparse matrix representation for large matrices with many zeros via --sparse flag
- Display results formatted as aligned matrices to console with configurable decimal precision via --precision flag (default: 4)
- Support both integer and floating-point matrix elements
- Verify results where possible: multiplying a matrix by its inverse should yield the identity matrix (within floating-point tolerance)
- Save results as JSON or CSV with --output flag (default: result_matrix.json)
- If no input is given, generate sample matrices and demonstrate each operation: add two 3x3 matrices, multiply a 3x3 by a 3x4, compute determinant and inverse of a 4x4, find eigenvalues of a symmetric 3x3, and show LU decomposition of a 4x4
- Handle errors: dimension mismatches, singular matrices (non-invertible), non-square matrices for operations requiring square input, and malformed input files

Create a complete C++ project for a clean Ubuntu 22.04 machine with only G++ 12+ and CMake 3.22+ installed. Include:
- Source code
- CMakeLists.txt with all dependencies (direct and transitive) pinned to exact versions
- README.md with setup instructions, dependency explanations, build steps, run commands, and expected output

[]

cmake_minimum_required(VERSION 3.16)
project(matrix_operations_calculator LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

add_executable(matrix_operations_calculator src/main.cpp)

# Matrix Operations Calculator (C++)

## Requirements
- Ubuntu 22.04
- CMake 3.16+
- C++20 compiler

## Build
```bash
cmake -S . -B build
cmake --build build
```

## Run
```bash
./build/matrix_operations_calculator matrix.csv --operation inverse --precision 6
./build/matrix_operations_calculator matrix.json --operation multiply --matrix2 matrix2.json
./build/matrix_operations_calculator matrix.csv --operation eigenvalues --output result_matrix.json
```

If no input is provided, sample operations are demonstrated.

#include <cmath>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>

using Matrix = std::vector<std::vector<double>>;
namespace fs = std::filesystem;

struct Config {
    std::string matrix1;
    std::string matrix2;
    std::string operation = "transpose";
    bool hasScalar = false;
    double scalar = 0.0;
    bool sparse = false;
    int precision = 4;
    std::string output = "result_matrix.json";
};

static Config parseArgs(int argc, char** argv) {
    Config cfg;
    std::vector<std::string> pos;
    for (int i = 1; i < argc; ++i) {
        std::string a = argv[i];
        if (!a.starts_with("--")) {
            pos.push_back(a);
            continue;
        }
        if (a == "--operation" && i + 1 < argc) cfg.operation = argv[++i];
        else if (a == "--matrix2" && i + 1 < argc) cfg.matrix2 = argv[++i];
        else if (a == "--scalar" && i + 1 < argc) { cfg.hasScalar = true; cfg.scalar = std::stod(argv[++i]); }
        else if (a == "--sparse") cfg.sparse = true;
        else if (a == "--precision" && i + 1 < argc) cfg.precision = std::stoi(argv[++i]);
        else if (a == "--output" && i + 1 < argc) cfg.output = argv[++i];
        else throw std::runtime_error("Unknown option: " + a);
    }
    if (!pos.empty()) cfg.matrix1 = pos[0];
    return cfg;
}

static Matrix readCsvMatrix(const fs::path& p) {
    std::ifstream in(p);
    if (!in) throw std::runtime_error("Cannot open file: " + p.string());
    Matrix m;
    std::string line;
    while (std::getline(in, line)) {
        if (line.empty()) continue;
        std::stringstream ss(line);
        std::string cell;
        std::vector<double> row;
        while (std::getline(ss, cell, ',')) row.push_back(std::stod(cell));
        m.push_back(row);
    }
    return m;
}

static Matrix readJsonMatrix(const fs::path& p) {
    std::string txt;
    {
        std::ifstream in(p);
        if (!in) throw std::runtime_error("Cannot open file: " + p.string());
        std::ostringstream ss;
        ss << in.rdbuf();
        txt = ss.str();
    }
    txt.erase(std::remove_if(txt.begin(), txt.end(), [](unsigned char c){ return std::isspace(c); }), txt.end());
    if (txt.size() < 4 || txt.substr(0, 2) != "[[" || txt.substr(txt.size() - 2) != "]]") throw std::runtime_error("Invalid JSON matrix");
    txt = txt.substr(2, txt.size() - 4);
    Matrix m;
    std::size_t start = 0;
    while (true) {
        std::size_t pos = txt.find("],[", start);
        std::string rowtxt = pos == std::string::npos ? txt.substr(start) : txt.substr(start, pos - start);
        std::stringstream ss(rowtxt);
        std::string v;
        std::vector<double> row;
        while (std::getline(ss, v, ',')) row.push_back(std::stod(v));
        m.push_back(row);
        if (pos == std::string::npos) break;
        start = pos + 3;
    }
    return m;
}

static Matrix readMatrix(const std::string& path) {
    fs::path p(path);
    if (p.extension() == ".json") return readJsonMatrix(p);
    return readCsvMatrix(p);
}

static std::string toJsonString(const std::string& s) {
    std::string out = "\"";
    for (char c : s) {
        if (c == '\\') out += "\\\\";
        else if (c == '"') out += "\\\"";
        else if (c == '\n') out += "\\n";
        else out += c;
    }
    out += "\"";
    return out;
}

static std::string matrixToJson(const Matrix& m) {
    std::ostringstream out;
    out << "[";
    for (std::size_t i = 0; i < m.size(); ++i) {
        if (i) out << ",";
        out << "[";
        for (std::size_t j = 0; j < m[i].size(); ++j) {
            if (j) out << ",";
            out << m[i][j];
        }
        out << "]";
    }
    out << "]";
    return out.str();
}

static std::string formatMatrix(const Matrix& m, int precision) {
    if (m.empty()) return "[]";
    std::vector<std::vector<std::string>> txt(m.size(), std::vector<std::string>(m[0].size()));
    std::vector<std::size_t> w(m[0].size(), 0);
    for (std::size_t i = 0; i < m.size(); ++i) {
        for (std::size_t j = 0; j < m[i].size(); ++j) {
            std::ostringstream s;
            s.setf(std::ios::fixed);
            s.precision(precision);
            s << m[i][j];
            txt[i][j] = s.str();
            w[j] = std::max(w[j], txt[i][j].size());
        }
    }
    std::ostringstream out;
    for (std::size_t i = 0; i < txt.size(); ++i) {
        if (i) out << "\n";
        for (std::size_t j = 0; j < txt[i].size(); ++j) {
            if (j) out << "  ";
            out << std::setw(static_cast<int>(w[j])) << txt[i][j];
        }
    }
    return out.str();
}

static bool sameDims(const Matrix& a, const Matrix& b) { return !a.empty() && !b.empty() && a.size() == b.size() && a[0].size() == b[0].size(); }
static bool canMul(const Matrix& a, const Matrix& b) { return !a.empty() && !b.empty() && a[0].size() == b.size(); }

static Matrix add(const Matrix& a, const Matrix& b, int sign = 1) {
    if (!sameDims(a, b)) throw std::runtime_error("Dimension mismatch for add/subtract");
    Matrix c = a;
    for (std::size_t i = 0; i < a.size(); ++i)
        for (std::size_t j = 0; j < a[0].size(); ++j)
            c[i][j] = a[i][j] + sign * b[i][j];
    return c;
}

static Matrix mul(const Matrix& a, const Matrix& b) {
    if (!canMul(a, b)) throw std::runtime_error("Dimension mismatch for multiply");
    Matrix c(a.size(), std::vector<double>(b[0].size(), 0.0));
    for (std::size_t i = 0; i < a.size(); ++i)
        for (std::size_t k = 0; k < a[0].size(); ++k)
            for (std::size_t j = 0; j < b[0].size(); ++j)
                c[i][j] += a[i][k] * b[k][j];
    return c;
}

static Matrix transpose(const Matrix& a) {
    Matrix t(a[0].size(), std::vector<double>(a.size()));
    for (std::size_t i = 0; i < a.size(); ++i)
        for (std::size_t j = 0; j < a[0].size(); ++j)
            t[j][i] = a[i][j];
    return t;
}

static double determinant(Matrix a) {
    std::size_t n = a.size();
    if (n == 0 || a[0].size() != n) throw std::runtime_error("Determinant requires square matrix");
    double det = 1.0;
    for (std::size_t i = 0; i < n; ++i) {
        std::size_t pivot = i;
        for (std::size_t r = i + 1; r < n; ++r) if (std::fabs(a[r][i]) > std::fabs(a[pivot][i])) pivot = r;
        if (std::fabs(a[pivot][i]) < 1e-12) return 0.0;
        if (pivot != i) {
            std::swap(a[pivot], a[i]);
            det = -det;
        }
        det *= a[i][i];
        for (std::size_t r = i + 1; r < n; ++r) {
            double factor = a[r][i] / a[i][i];
            for (std::size_t c = i; c < n; ++c) a[r][c] -= factor * a[i][c];
        }
    }
    return det;
}

static Matrix inverse(Matrix a) {
    std::size_t n = a.size();
    if (n == 0 || a[0].size() != n) throw std::runtime_error("Inverse requires square matrix");
    Matrix inv(n, std::vector<double>(n, 0.0));
    for (std::size_t i = 0; i < n; ++i) inv[i][i] = 1.0;

    for (std::size_t i = 0; i < n; ++i) {
        std::size_t pivot = i;
        for (std::size_t r = i + 1; r < n; ++r) if (std::fabs(a[r][i]) > std::fabs(a[pivot][i])) pivot = r;
        if (std::fabs(a[pivot][i]) < 1e-12) throw std::runtime_error("Singular matrix");
        if (pivot != i) {
            std::swap(a[pivot], a[i]);
            std::swap(inv[pivot], inv[i]);
        }
        double diag = a[i][i];
        for (std::size_t c = 0; c < n; ++c) { a[i][c] /= diag; inv[i][c] /= diag; }

        for (std::size_t r = 0; r < n; ++r) {
            if (r == i) continue;
            double factor = a[r][i];
            for (std::size_t c = 0; c < n; ++c) {
                a[r][c] -= factor * a[i][c];
                inv[r][c] -= factor * inv[i][c];
            }
        }
    }
    return inv;
}

static int rank(Matrix a) {
    std::size_t m = a.size(), n = a[0].size();
    std::size_t r = 0;
    for (std::size_t c = 0; c < n && r < m; ++c) {
        std::size_t pivot = r;
        for (std::size_t i = r + 1; i < m; ++i) if (std::fabs(a[i][c]) > std::fabs(a[pivot][c])) pivot = i;
        if (std::fabs(a[pivot][c]) < 1e-10) continue;
        std::swap(a[pivot], a[r]);
        double div = a[r][c];
        for (std::size_t j = c; j < n; ++j) a[r][j] /= div;
        for (std::size_t i = 0; i < m; ++i) {
            if (i == r) continue;
            double f = a[i][c];
            for (std::size_t j = c; j < n; ++j) a[i][j] -= f * a[r][j];
        }
        r++;
    }
    return static_cast<int>(r);
}

static double trace(const Matrix& a) {
    if (a.size() != a[0].size()) throw std::runtime_error("Trace requires square matrix");
    double t = 0;
    for (std::size_t i = 0; i < a.size(); ++i) t += a[i][i];
    return t;
}

static std::pair<Matrix, Matrix> lu(const Matrix& a) {
    std::size_t n = a.size();
    if (n == 0 || a[0].size() != n) throw std::runtime_error("LU requires square matrix");
    Matrix L(n, std::vector<double>(n, 0.0));
    Matrix U(n, std::vector<double>(n, 0.0));

    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t k = i; k < n; ++k) {
            double sum = 0;
            for (std::size_t j = 0; j < i; ++j) sum += L[i][j] * U[j][k];
            U[i][k] = a[i][k] - sum;
        }
        for (std::size_t k = i; k < n; ++k) {
            if (i == k) L[i][i] = 1;
            else {
                double sum = 0;
                for (std::size_t j = 0; j < i; ++j) sum += L[k][j] * U[j][i];
                if (std::fabs(U[i][i]) < 1e-12) throw std::runtime_error("LU failed due to zero pivot");
                L[k][i] = (a[k][i] - sum) / U[i][i];
            }
        }
    }
    return {L, U};
}

static bool isSymmetric(const Matrix& a) {
    if (a.size() != a[0].size()) return false;
    for (std::size_t i = 0; i < a.size(); ++i)
        for (std::size_t j = i + 1; j < a.size(); ++j)
            if (std::fabs(a[i][j] - a[j][i]) > 1e-9) return false;
    return true;
}

static bool isPositiveDefinite(const Matrix& a) {
    if (!isSymmetric(a)) return false;
    std::size_t n = a.size();
    Matrix L(n, std::vector<double>(n, 0.0));
    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = 0; j <= i; ++j) {
            double sum = 0;
            for (std::size_t k = 0; k < j; ++k) sum += L[i][k] * L[j][k];
            if (i == j) {
                double v = a[i][i] - sum;
                if (v <= 0) return false;
                L[i][j] = std::sqrt(v);
            } else {
                L[i][j] = (a[i][j] - sum) / L[j][j];
            }
        }
    }
    return true;
}

static double froNorm(const Matrix& a) {
    double s = 0;
    for (auto& r : a) for (double v : r) s += v * v;
    return std::sqrt(s);
}

static std::vector<double> eigenvaluesQR(Matrix a, int iterations = 200) {
    if (a.size() != a[0].size()) throw std::runtime_error("Eigenvalues require square matrix");
    std::size_t n = a.size();

    for (int it = 0; it < iterations; ++it) {
        // Classical Gram-Schmidt QR
        Matrix Q(n, std::vector<double>(n, 0.0));
        Matrix R(n, std::vector<double>(n, 0.0));
        Matrix V = a;
        for (std::size_t j = 0; j < n; ++j) {
            std::vector<double> v(n);
            for (std::size_t i = 0; i < n; ++i) v[i] = V[i][j];
            for (std::size_t k = 0; k < j; ++k) {
                double dot = 0;
                for (std::size_t i = 0; i < n; ++i) dot += Q[i][k] * v[i];
                R[k][j] = dot;
                for (std::size_t i = 0; i < n; ++i) v[i] -= dot * Q[i][k];
            }
            double norm = 0;
            for (double x : v) norm += x * x;
            norm = std::sqrt(norm);
            if (norm < 1e-12) norm = 1e-12;
            R[j][j] = norm;
            for (std::size_t i = 0; i < n; ++i) Q[i][j] = v[i] / norm;
        }
        a = mul(R, Q);
    }
    std::vector<double> vals(n);
    for (std::size_t i = 0; i < n; ++i) vals[i] = a[i][i];
    return vals;
}

static Matrix scalarMul(const Matrix& a, double s) {
    Matrix r = a;
    for (auto& row : r) for (double& v : row) v *= s;
    return r;
}

int main(int argc, char** argv) {
    try {
        Config cfg = parseArgs(argc, argv);
        if (cfg.matrix1.empty()) {
            Matrix A = {{4,2,1},{0,5,3},{2,1,6}};
            Matrix B = {{1,2,3},{4,5,6},{7,8,9}};
            auto [L,U] = lu(A);
            std::ostringstream demo;
            demo << "{\"add\":" << matrixToJson(add(A,B))
                 << ",\"determinant\":" << determinant(A)
                 << ",\"inverse\":" << matrixToJson(inverse(A))
                 << ",\"eigenvalues\":[";
            auto eig = eigenvaluesQR(A);
            for (std::size_t i=0;i<eig.size();++i){ if(i)demo<<","; demo<<eig[i]; }
            demo << "],\"lu\":{\"L\":" << matrixToJson(L) << ",\"U\":" << matrixToJson(U) << "}}";
            std::ofstream(cfg.output) << demo.str();
            std::cout << demo.str() << "\n";
            return 0;
        }

        Matrix A = readMatrix(cfg.matrix1);
        Matrix B;
        if (!cfg.matrix2.empty()) B = readMatrix(cfg.matrix2);

        std::ostringstream report;
        report << "{\n"
               << "  \"operation\": " << toJsonString(cfg.operation) << ",\n"
               << "  \"properties\": {\n"
               << "    \"dimensions\": [" << A.size() << "," << (A.empty()?0:A[0].size()) << "],\n"
               << "    \"rank\": " << rank(A) << ",\n"
               << "    \"trace\": " << ((A.size()==A[0].size()) ? std::to_string(trace(A)) : "null") << ",\n"
               << "    \"isSymmetric\": " << (isSymmetric(A)?"true":"false") << ",\n"
               << "    \"isPositiveDefinite\": " << (isPositiveDefinite(A)?"true":"false") << ",\n";
        if (A.size()==A[0].size()) {
            Matrix invA = inverse(A);
            report << "    \"conditionNumber\": " << (froNorm(A)*froNorm(invA)) << "\n";
        } else {
            report << "    \"conditionNumber\": null\n";
        }
        report << "  },\n";

        if (cfg.hasScalar) report << "  \"scalarMultiply\": " << matrixToJson(scalarMul(A,cfg.scalar)) << ",\n";

        if (cfg.operation == "add") report << "  \"result\": " << matrixToJson(add(A,B)) << "\n";
        else if (cfg.operation == "subtract") report << "  \"result\": " << matrixToJson(add(A,B,-1)) << "\n";
        else if (cfg.operation == "multiply") report << "  \"result\": " << matrixToJson(mul(A,B)) << "\n";
        else if (cfg.operation == "transpose") report << "  \"result\": " << matrixToJson(transpose(A)) << "\n";
        else if (cfg.operation == "determinant") report << "  \"result\": " << determinant(A) << "\n";
        else if (cfg.operation == "inverse") {
            Matrix invA = inverse(A);
            report << "  \"result\": " << matrixToJson(invA) << ",\n";
            report << "  \"inverseVerification\": " << matrixToJson(mul(A, invA)) << "\n";
        }
        else if (cfg.operation == "eigenvalues") {
            auto vals = eigenvaluesQR(A);
            report << "  \"result\": [";
            for (std::size_t i = 0; i < vals.size(); ++i) { if (i) report << ","; report << vals[i]; }
            report << "]\n";
        }
        else if (cfg.operation == "rank") report << "  \"result\": " << rank(A) << "\n";
        else if (cfg.operation == "trace") report << "  \"result\": " << trace(A) << "\n";
        else if (cfg.operation == "lu") {
            auto [L,U] = lu(A);
            report << "  \"result\": {\"L\": " << matrixToJson(L) << ", \"U\": " << matrixToJson(U) << "}\n";
        }
        else throw std::runtime_error("Unsupported operation: " + cfg.operation);

        report << "}\n";

        std::ofstream(cfg.output) << report.str();

        // Console display
        if (cfg.operation == "add") std::cout << formatMatrix(add(A,B), cfg.precision) << "\n";
        else if (cfg.operation == "subtract") std::cout << formatMatrix(add(A,B,-1), cfg.precision) << "\n";
        else if (cfg.operation == "multiply") std::cout << formatMatrix(mul(A,B), cfg.precision) << "\n";
        else if (cfg.operation == "transpose") std::cout << formatMatrix(transpose(A), cfg.precision) << "\n";
        else if (cfg.operation == "inverse") std::cout << formatMatrix(inverse(A), cfg.precision) << "\n";
        else std::cout << report.str() << "\n";
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << "\n";
        return 1;
    }
    return 0;
}