WIP

algebra/builtin/builtin_linalg.c

@@ -8,11 +8,12 @@
  * This source code is licensed under the BSD 3-Clause License
  */
 
-#include "qoco_linalg.h"
+#include "builtin_types.h"
 
-QOCOCscMatrix* new_qoco_csc_matrix(const QOCOCscMatrix* A)
+QOCOMatrix* new_qoco_matrix(const QOCOCscMatrix* A)
 {
-  QOCOCscMatrix* M = qoco_malloc(sizeof(QOCOCscMatrix));
+  QOCOMatrix* M = qoco_malloc(sizeof(QOCOMatrix));
+  QOCOCscMatrix* Mcsc = qoco_malloc(sizeof(QOCOCscMatrix));
 
   if (A) {
     QOCOInt m = A->m;
@@ -27,61 +28,88 @@ QOCOCscMatrix* new_qoco_csc_matrix(const QOCOCscMatrix* A)
     copy_arrayi(A->i, i, nnz);
     copy_arrayi(A->p, p, n + 1);
 
-    M->m = m;
-    M->n = n;
-    M->nnz = nnz;
-    M->x = x;
-    M->i = i;
-    M->p = p;
+    Mcsc->m = m;
+    Mcsc->n = n;
+    Mcsc->nnz = nnz;
+    Mcsc->x = x;
+    Mcsc->i = i;
+    Mcsc->p = p;
   }
   else {
-    M->m = 0;
-    M->n = 0;
-    M->nnz = 0;
-    M->x = NULL;
-    M->i = NULL;
-    M->p = NULL;
+    Mcsc->m = 0;
+    Mcsc->n = 0;
+    Mcsc->nnz = 0;
+    Mcsc->x = NULL;
+    Mcsc->i = NULL;
+    Mcsc->p = NULL;
   }
 
+  M->csc = M;
+
   return M;
 }
 
-void free_qoco_csc_matrix(QOCOCscMatrix* A)
+QOCOVectorf* new_qoco_vectorf(const QOCOFloat* x, QOCOInt n)
 {
-  free(A->x);
-  free(A->i);
-  free(A->p);
-  free(A);
+  QOCOVectorf* v = qoco_malloc(sizeof(QOCOVectorf));
+  QOCOFloat* vdata = qoco_malloc(sizeof(QOCOFloat) * n);
+  copy_arrayf(x, vdata, n);
+
+  v->len = n;
+  v->data = vdata;
 }
 
-void copy_arrayf(const QOCOFloat* x, QOCOFloat* y, QOCOInt n)
+void free_qoco_matrix(QOCOMatrix* A)
 {
-  qoco_assert(x || n == 0);
-  qoco_assert(y || n == 0);
+  free_qoco_csc_matrix(A);
+  qoco_free(A);
+}
 
-  for (QOCOInt i = 0; i < n; ++i) {
-    y[i] = x[i];
-  }
+void free_qoco_vectorf(QOCOVectorf* x)
+{
+  qoco_free(x->data);
+  qoco_free(x);
 }
 
-void copy_and_negate_arrayf(const QOCOFloat* x, QOCOFloat* y, QOCOInt n)
+QOCOMatrix* create_transposed_matrix(const QOCOMatrix* A, QOCOInt* AtoAt)
 {
-  qoco_assert(x || n == 0);
-  qoco_assert(y || n == 0);
+  QOCOMatrix* At = qoco_malloc(sizeof(QOCOMatrix));
+  QOCOCscMatrix* Atcsc = create_transposed_csc_matrix(A->csc, AtoAt);
+  At->csc = Atcsc;
+  return At;
+}
 
-  for (QOCOInt i = 0; i < n; ++i) {
-    y[i] = -x[i];
-  }
+QOCOMatrix* regularize_P(QOCOInt num_diagP, QOCOMatrix* P, QOCOFloat reg,
+                         QOCOInt* nzadded_idx)
+{
+  regularize_P_csc(num_diagP, P->csc, reg, nzadded_idx);
+  qoco_free(P);
 }
 
-void copy_arrayi(const QOCOInt* x, QOCOInt* y, QOCOInt n)
+void unregularize(QOCOMatrix* M, QOCOFloat lambda)
 {
-  qoco_assert(x || n == 0);
-  qoco_assert(y || n == 0);
+  unregularize_csc(M->csc, lambda);
+}
 
-  for (QOCOInt i = 0; i < n; ++i) {
-    y[i] = x[i];
-  }
+QOCOMatrix* construct_identity(QOCOInt n, QOCOFloat lambda)
+{
+  QOCOMatrix* M = qoco_malloc(sizeof(QOCOMatrix));
+  M->csc = construct_identity_csc(n, lambda);
+}
+
+void scale_matrix(QOCOFloat a, QOCOMatrix* M)
+{
+  scale_arrayf(M->csc->x, M->csc->x, a, M->csc->nnz);
+}
+
+void row_col_scale(QOCOMatrix* M, QOCOVectorf* E, QOCOVectorf* D)
+{
+  row_col_scale_csc(M->csc, E->data, D->data);
+}
+
+void update_matrix(QOCOMatrix* M, QOCOFloat* Mnew)
+{
+  copy_arrayf(Mnew, M->csc->x, M->csc->nnz);
 }
 
 QOCOFloat qoco_dot(const QOCOFloat* u, const QOCOFloat* v, QOCOInt n)
@@ -96,24 +124,33 @@ QOCOFloat qoco_dot(const QOCOFloat* u, const QOCOFloat* v, QOCOInt n)
   return x;
 }
 
-QOCOInt max_arrayi(const QOCOInt* x, QOCOInt n)
+void ew_product(QOCOVectorf* x, QOCOVectorf* y, QOCOVectorf* z)
 {
-  qoco_assert(x || n == 0);
+  qoco_assert(x->len == y->len);
+  qoco_assert(x->len == z->len);
+  ew_product_arrayf(x->data, y->data, z->data, x->len);
+}
 
-  QOCOInt max = -QOCOInt_MAX;
-  for (QOCOInt i = 0; i < n; ++i) {
-    max = qoco_max(max, x[i]);
-  }
-  return max;
+void ew_product_vec_array(QOCOVectorf* x, QOCOFloat* y, QOCOVectorf* z)
+{
+  qoco_assert(x->len == y->len);
+  qoco_assert(x->len == z->len);
+  ew_product_arrayf(x->data, y, z->data, x->len);
 }
 
-void scale_arrayf(const QOCOFloat* x, QOCOFloat* y, QOCOFloat s, QOCOInt n)
+void scale_vectorf(QOCOFloat a, QOCOVectorf* u)
 {
-  qoco_assert(x || n == 0);
-  qoco_assert(y || n == 0);
+  scale_arrayf(u->data, u->data, a, u->len);
+}
 
-  for (QOCOInt i = 0; i < n; ++i) {
-    y[i] = s * x[i];
+void copy_vectorf(QOCOVectorf* src, QOCOFloat* dest, QOCOInt dest_idx,
+                  QOCOInt negate)
+{
+  if (negate) {
+    copy_arrayf(src->data, &dest[dest_idx], src->len);
+  }
+  else {
+    copy_and_negate_arrayf(src->data, &dest[dest_idx], src->len);
   }
 }
 

algebra/builtin/builtin_types.h

@@ -0,0 +1,36 @@
+/**
+ * @file builtin_types.h
+ * @author Govind M. Chari <govindchari1@gmail.com>
+ *
+ * @section LICENSE
+ *
+ * Copyright (c) 2025, Govind M. Chari
+ * This source code is licensed under the BSD 3-Clause License
+ *
+ * @section DESCRIPTION
+ *
+ * Defines the vector and matrices for builtin linear algebra.
+ */
+
+#ifndef BUILTIN_TYPES_H
+#define BUILTIN_TYPES_H
+
+#include "definitions.h"
+#include "qoco_linalg.h"
+#include "common_linalg.h"
+
+struct QOCOVectori_ {
+  QOCOInt* data;
+  QOCOInt  len;
+};
+
+struct QOCOVectorf_ {
+  QOCOFloat* data;
+  QOCOInt    len;
+};
+
+struct QOCOMatrix_ {
+  QOCOCscMatrix* csc;
+};
+
+#endif /* ifndef BUILTIN_TYPES_H */

devtools/run_tests.sh

@@ -1 +1 @@
-export CXX=/usr/local/bin/clang++ && export CC=/usr/local/bin/clang && cd build && cmake -DQOCO_BUILD_TYPE:STR=Release -DENABLE_TESTING:BOOL=True -DBUILD_QOCO_BENCHMARK_RUNNER:BOOL=True .. && make && ctest --verbose && cd ..
+export CXX=/usr/bin/clang++ && export CC=/usr/bin/clang && cd build && cmake -DQOCO_BUILD_TYPE:STR=Release -DENABLE_TESTING:BOOL=True -DBUILD_QOCO_BENCHMARK_RUNNER:BOOL=True .. && make && ctest --verbose && cd ..

include/common_linalg.h

@@ -17,14 +17,22 @@
 #include "definitions.h"
 #include "qoco_linalg.h"
 
+/**
+ * @brief Frees all the internal arrays and the pointer to the QOCOCscMatrix.
+ * Should only be used if QOCOCscMatrix and all internal arrays were malloc'ed.
+ *
+ * @param A Pointer to QOCOCscMatrix.
+ */
+void free_qoco_csc_matrix(QOCOCscMatrix* A);
+
 /**
  * @brief Allocates a new csc matrix that is lambda * I.
  *
  * @param n Size of identity matrix.
  * @param lambda Scaling factor for identity.
  * @return Pointer to new constructed matrix.
  */
-QOCOCscMatrix* construct_identity(QOCOInt n, QOCOFloat lambda);
+QOCOCscMatrix* construct_identity_csc(QOCOInt n, QOCOFloat lambda);
 
 /**
  * @brief Counts the number of diagonal elements in upper triangular CSC matrix
@@ -46,7 +54,7 @@ QOCOInt count_diag(QOCOCscMatrix* M);
  * @param nzadded_idx Indices of elements of M->x that are added.
  * @return P + reg * I.
  */
-QOCOCscMatrix* regularize_P(QOCOInt num_diagP, QOCOCscMatrix* P, QOCOFloat reg,
+QOCOCscMatrix* regularize_P_csc(QOCOInt num_diagP, QOCOCscMatrix* P, QOCOFloat reg,
                             QOCOInt* nzadded_idx);
 
 /**
@@ -57,7 +65,7 @@ QOCOCscMatrix* regularize_P(QOCOInt num_diagP, QOCOCscMatrix* P, QOCOFloat reg,
  * @param M Matrix.
  * @param lambda Regularization.
  */
-void unregularize(QOCOCscMatrix* M, QOCOFloat lambda);
+void unregularize_csc(QOCOCscMatrix* M, QOCOFloat lambda);
 
 /**
  * @brief Computes the infinity norm of each column (or equivalently row) of a
@@ -83,17 +91,44 @@ void row_inf_norm(const QOCOCscMatrix* M, QOCOFloat* norm);
  * @param A Input matrix.
  * @param AtoAt Mapping from A to At.
  */
-QOCOCscMatrix* create_transposed_matrix(const QOCOCscMatrix* A, QOCOInt* AtoAt);
+QOCOCscMatrix* create_transposed_csc_matrix(const QOCOCscMatrix* A, QOCOInt* AtoAt);
 
 /**
  * @brief Scales the rows of M by E and columns of M by D.
  * M = diag(E) * M * diag(S)
  *
- * @param M An m by n sparse matrix.
+ * @param M An m by n QOCOCscMatrix matrix.
  * @param E Vector of length m.
  * @param D Vector of length m.
  */
-void row_col_scale(const QOCOCscMatrix* M, QOCOFloat* E, QOCOFloat* D);
+void row_col_scale_csc(const QOCOCscMatrix* M, QOCOFloat* E, QOCOFloat* D);
+
+/**
+ * @brief Copies array of QOCOFloats from x to array y.
+ *
+ * @param x Source array.
+ * @param y Destination array.
+ * @param n Length of arrays.
+ */
+void copy_arrayf(const QOCOFloat* x, QOCOFloat* y, QOCOInt n);
+
+/**
+ * @brief Copies and negates array of QOCOFloats from x to array y.
+ *
+ * @param x Source array.
+ * @param y Destination array.
+ * @param n Length of arrays.
+ */
+void copy_and_negate_arrayf(const QOCOFloat* x, QOCOFloat* y, QOCOInt n);
+
+/**
+ * @brief Copies array of QOCOInts from x to array y.
+ *
+ * @param x Source array.
+ * @param y Destination array.
+ * @param n Length of arrays.
+ */
+void copy_arrayi(const QOCOInt* x, QOCOInt* y, QOCOInt n);
 
 /**
  * @brief Computes elementwise product z = x .* y
@@ -103,7 +138,27 @@ void row_col_scale(const QOCOCscMatrix* M, QOCOFloat* E, QOCOFloat* D);
  * @param z Output array.
  * @param n Length of arrays.
  */
-void ew_product(QOCOFloat* x, const QOCOFloat* y, QOCOFloat* z, QOCOInt n);
+void ew_product_arrayf(QOCOFloat* x, QOCOFloat* y, QOCOFloat* z, QOCOInt n);
+
+/**
+ * @brief Computes maximum element of array of QOCOInts.
+ *
+ * @param x Input array.
+ * @param n Length of array.
+ * @return Maximum element of x.
+ */
+QOCOInt max_arrayi(const QOCOInt* x, QOCOInt n);
+
+/**
+ * @brief Scales array x by s and stores result in y.
+ * y = s * x
+ *
+ * @param x Input array.
+ * @param y Output array.
+ * @param s Scaling factor.
+ * @param n Length of arrays.
+ */
+void scale_arrayf(QOCOFloat* x, QOCOFloat* y, QOCOFloat s, QOCOInt n);
 
 /**
  * @brief Inverts permutation vector p and stores inverse in pinv.

include/kkt.h

@@ -21,6 +21,12 @@
 #include "qoco_linalg.h"
 #include "structs.h"
 
+QOCOCscMatrix* create_kkt(QOCOCscMatrix* P, QOCOCscMatrix* A, QOCOCscMatrix* G,
+                          QOCOCscMatrix* At, QOCOCscMatrix* Gt, QOCOFloat static_reg, QOCOInt n,
+                          QOCOInt m, QOCOInt p, QOCOInt l, QOCOInt nsoc,
+                          QOCOInt* q, QOCOInt* PregtoKKT, QOCOInt* AttoKKT,
+                          QOCOInt* GttoKKT, QOCOInt* nt2kkt, QOCOInt* ntdiag2kkt);
+
 /**
  * @brief Allocate memory for KKT matrix.
  *