From 880fef6d22f2d2382765b10a843355afebdd5335 Mon Sep 17 00:00:00 2001 From: araujoms Date: Wed, 25 Dec 2024 11:10:33 +0100 Subject: [PATCH] add complex psd cone --- docs/src/api/cones.rst | 65 ++++++---- include/cones.h | 6 +- include/scs.h | 4 + include/scs_blas.h | 4 + include/scs_types.h | 4 + src/cones.c | 268 +++++++++++++++++++++++++++++++++++++++-- 6 files changed, 321 insertions(+), 30 deletions(-) diff --git a/docs/src/api/cones.rst b/docs/src/api/cones.rst index aced3b32..522a4784 100644 --- a/docs/src/api/cones.rst +++ b/docs/src/api/cones.rst @@ -49,6 +49,10 @@ The cone :math:`\mathcal{K}` can be any Cartesian product of the following primi - :math:`\{ (u,v,w)\in \mathbf{R}^3 \mid \left(\frac{u}{p}\right)^p \left(\frac{v}{1-p}\right)^{1-p} \geq |w|\}` - :code:`p` array :math:`p\in[-1,1]` powers with :code:`psize` elements, negative entries correspond to dual power cone (sign is flipped). - :math:`3 \times`:code:`psize` (total for primal and dual power cone) + * - Complex positive semidefinite cone + - :math:`\{ s \in \mathbf{R}^{k^2} \mid \text{mat}(s) \succeq 0 \}` (See :ref:`note `) + - :code:`cs` array of PSD cone lengths with :code:`cssize` elements, each :math:`cs_i = k_i` in description. + - :math:`\displaystyle \sum_{i=1}^{\text{cssize}} cs_i^2` **Note**: @@ -70,27 +74,30 @@ and therefore the variable vector is stacked as :code:`[x, y, z]`, etc. Semidefinite cones ------------------ -The symmetric positive semidefinite cone of matrices is the set +The real positive semidefinite cone of matrices is the set .. math:: - \{S \in \mathbf{R}^{k \times k} \mid S = S^\top, x^\top S x \geq 0 \ \forall x \in \mathbf{R}^k \} + \{S \in \mathbf{R}^{k \times k} \mid S = S^\top,\ x^\top S x \geq 0 \ \forall x \in \mathbf{R}^k \}, + +and the complex positive semidefinite cone is the set + +.. math:: + \{S \in \mathbf{C}^{k \times k} \mid S = S^\dagger,\ x^\dagger S x \geq 0 \ \forall x \in \mathbf{C}^k \}, and for short we use :math:`S \succeq 0` to denote membership. SCS -vectorizes this cone in a special way which we detail here. +vectorizes these cones in a special way which we detail here. -SCS assumes that the input data corresponding to -semidefinite cones have been vectorized by scaling the off-diagonal entries by -:math:`\sqrt{2}` and stacking the lower triangular elements column-wise. For a :math:`k \times k` -matrix variable (or data matrix) this operation would create a vector of length -:math:`k(k+1)/2`. Scaling by :math:`\sqrt{2}` is required to preserve the inner-product. +SCS assumes that the input data corresponding to positive semidefinite cones have been vectorized by scaling the off-diagonal entries by +:math:`\sqrt{2}` and stacking the real parameters of the upper triangular row-wise. For a :math:`k \times k` +real matrix variable (or data matrix) this operation would create a vector of length +:math:`k(k+1)/2`, whereas for a :math:`k \times k` complex matrix its creates a vector of length :math:`k^2`. Scaling by :math:`\sqrt{2}` is required to preserve the inner-product. **This must be done for the rows of both** :math:`A` **and** :math:`b` **that correspond to semidefinite cones and must be done independently for each semidefinite cone.** -More explicitly, we want to express :math:`\text{Trace}(Y S)` as :math:`\text{vec}(Y)^\top \text{vec}(S)`, -where the :math:`\text{vec}` operation takes the (assumed to be symmetric) :math:`k \times k` matrix +More explicitly, we want to express :math:`\text{Trace}(Y^\dagger S)` as :math:`\text{vec}(Y)^\dagger \text{vec}(S)`, +where the :math:`\text{vec}` operation takes the :math:`k \times k` matrix .. math:: - S = \begin{bmatrix} S_{11} & S_{12} & \ldots & S_{1k} \\ S_{21} & S_{22} & \ldots & S_{2k} \\ @@ -98,17 +105,21 @@ where the :math:`\text{vec}` operation takes the (assumed to be symmetric) :math S_{k1} & S_{k2} & \ldots & S_{kk} \\ \end{bmatrix} -and produces a vector consisting of the lower triangular elements scaled and arranged as +and produces a vector consisting of the upper triangular real parameters elements scaled and re-arranged. In the real case, the result is .. math:: + \text{vec}(S) = (S_{11}, \sqrt{2} S_{12}, \ldots, \sqrt{2} S_{1k}, S_{22}, \sqrt{2}S_{23}, \dots, S_{k-1,k-1}, \sqrt{2}S_{k-1,k}, S_{kk}) \in \mathbf{R}^{k(k+1)/2}, + +whereas in the complex case the result is - \text{vec}(S) = (S_{11}, \sqrt{2} S_{21}, \ldots, \sqrt{2} S_{k1}, S_{22}, \sqrt{2}S_{32}, \dots, S_{k-1,k-1}, \sqrt{2}S_{k,k-1}, S_{kk}) \in \mathbf{R}^{k(k+1)/2}. +.. math:: + \text{vec}(S) = (S_{11}, \sqrt{2} \Re(S_{12}), \sqrt{2} \Im(S_{12}), \ldots, \sqrt{2} \Re(S_{1k}), \sqrt{2} \Im(S_{1k}), S_{22}, \\\sqrt{2}\Re(S_{23}), \sqrt{2}\Im(S_{23}), \dots, S_{k-1,k-1}, \sqrt{2}\Re(S_{k-1,k}), \sqrt{2}\Im(S_{k-1,k}), S_{kk}) \in \mathbf{R}^{k^2}. To recover the matrix solution this operation must be inverted on the components of the vectors returned by SCS corresponding to each semidefinite cone. That is, the -off-diagonal entries must be scaled by :math:`1/\sqrt{2}` and the upper triangular -entries are filled in by copying the values of lower triangular entries. -Explicitly, the inverse operation takes vector :math:`s \in +off-diagonal entries must be scaled by :math:`1/\sqrt{2}` and the lower triangular +entries are filled in by copying the values of upper triangular entries. +Explicitly, in the real case the inverse operation takes vector :math:`s \in \mathbf{R}^{k(k+1)/2}` and produces the matrix .. math:: @@ -118,19 +129,31 @@ Explicitly, the inverse operation takes vector :math:`s \in \vdots & \vdots & \ddots & \vdots \\ s_{k} / \sqrt{2} & s_{2k-1} / \sqrt{2} & \ldots & s_{k(k+1) / 2} \\ \end{bmatrix} - \in \mathbf{R}^{k \times k}. + \in \mathbf{R}^{k \times k}, +whereas in the complex case the inverse operation takes vector :math:`s \in +\mathbf{R}^{k^2}` and produces the matrix + +.. math:: + \text{mat}(s) = \begin{bmatrix} + s_{1} & (s_{2} + i s_3) / \sqrt{2} & \ldots & (s_{2k-2}+is_{2k-1}) / \sqrt{2} \\ + (s_{2} - i s_3) / \sqrt{2} / \sqrt{2} & s_{2k} & \ldots & (s_{4k-5}+is_{4k-4}) / \sqrt{2} \\ + \vdots & \vdots & \ddots & \vdots \\ + (s_{2k-2}-is_{2k-1}) / \sqrt{2} & (s_{4k-5}-is_{4k-4}) / \sqrt{2} & \ldots & s_{k^2} \\ + \end{bmatrix} + \in \mathbf{C}^{k \times k}, -So the cone definition that SCS uses is +So the cone definitions that SCS uses are .. math:: - \mathcal{S}_+^k = \{ \text{vec}(S) \mid S \succeq 0\} = \{s \in \mathbf{R}^{k(k+1)/2} \mid \text{mat}(s) \succeq 0 \}. + \mathcal{S}_\mathbf{R}^k = \{ \text{vec}(S) \mid S \succeq 0\} = \{s \in \mathbf{R}^{k(k+1)/2} \mid \text{mat}(s) \succeq 0 \}.\\ + \mathcal{S}_\mathbf{C}^k = \{ \text{vec}(S) \mid S \succeq 0\} = \{s \in \mathbf{R}^{k^2 } \mid \text{mat}(s) \succeq 0 \}. Example ^^^^^^^ -For a concrete example in python see :ref:`py_mat_completion`. -Here we consider the symmetric positive semidefinite cone constraint over +For a concrete example in Python see :ref:`py_mat_completion`. +Here we consider the real positive semidefinite cone constraint over variables :math:`x \in \mathbf{R}^n` and :math:`S \in \mathbf{R}^{k \times k}` .. math:: diff --git a/include/cones.h b/include/cones.h index 44c29611..8557f709 100644 --- a/include/cones.h +++ b/include/cones.h @@ -28,8 +28,10 @@ struct SCS_CONE_WORK { scs_float box_t_warm_start; #ifdef USE_LAPACK /* workspace for eigenvector decompositions: */ - scs_float *Xs, *Z, *e, *work; - blas_int lwork; + scs_float *Xs, *Z, *e, *work, *rwork; + scs_complex_float *cXs, *cZ, *cwork; + blas_int *isuppz, *iwork; + blas_int lwork, lcwork, lrwork, liwork; #endif }; diff --git a/include/scs.h b/include/scs.h index 4059a018..7042c77a 100644 --- a/include/scs.h +++ b/include/scs.h @@ -130,6 +130,10 @@ typedef struct { scs_int *s; /** Length of semidefinite constraints array `s`. */ scs_int ssize; + /** Array of complex semidefinite cone constraints, `len(cs) = cssize`. */ + scs_int *cs; + /** Length of complex semidefinite constraints array `cs`. */ + scs_int cssize; /** Number of primal exponential cone triples. */ scs_int ep; /** Number of dual exponential cone triples. */ diff --git a/include/scs_blas.h b/include/scs_blas.h index ed9807d7..15488904 100644 --- a/include/scs_blas.h +++ b/include/scs_blas.h @@ -18,9 +18,11 @@ extern "C" { #ifndef SFLOAT #define BLAS(x) d##x #define BLASI(x) id##x +#define BLASC(x) z##x #else #define BLAS(x) s##x #define BLASI(x) is##x +#define BLASC(x) c##x #endif #else /* this extra indirection is needed for BLASSUFFIX to work correctly as a @@ -31,9 +33,11 @@ extern "C" { #ifndef SFLOAT #define BLAS(x) stitch__(d, x, BLASSUFFIX) #define BLASI(x) stitch__(id, x, BLASSUFFIX) +#define BLASC(x) stitch__(z, x, BLASSUFFIX) #else #define BLAS(x) stitch__(s, x, BLASSUFFIX) #define BLASI(x) stitch__(is, x, BLASSUFFIX) +#define BLASC(x) stitch__(c, x, BLASSUFFIX) #endif #endif diff --git a/include/scs_types.h b/include/scs_types.h index 02cfc27d..f65a6a68 100644 --- a/include/scs_types.h +++ b/include/scs_types.h @@ -11,6 +11,8 @@ extern "C" { #endif +#include + #ifdef DLONG /*#ifdef _WIN64 #include @@ -26,8 +28,10 @@ typedef int scs_int; #ifndef SFLOAT typedef double scs_float; +typedef double complex scs_complex_float; #else typedef float scs_float; +typedef float complex scs_complex_float; #endif #ifdef __cplusplus diff --git a/src/cones.c b/src/cones.c index 0d50a7ce..e8b7eaef 100644 --- a/src/cones.c +++ b/src/cones.c @@ -20,12 +20,24 @@ extern "C" { void BLAS(syev)(const char *jobz, const char *uplo, blas_int *n, scs_float *a, blas_int *lda, scs_float *w, scs_float *work, blas_int *lwork, blas_int *info); +void BLASC(heevr)(const char *jobz, const char *range, const char *uplo, blas_int *n, scs_complex_float *a, + blas_int *lda, scs_float *vl, scs_float *vu, blas_int *il, blas_int *iu, scs_float *abstol, + blas_int *m, scs_float *w, scs_complex_float *z, blas_int *ldz, blas_int *isuppz, scs_complex_float *cwork, + blas_int *lcwork, scs_float *rwork, blas_int *lrwork, blas_int *iwork, blas_int *liwork, blas_int *info); + blas_int BLAS(syrk)(const char *uplo, const char *trans, const blas_int *n, const blas_int *k, const scs_float *alpha, const scs_float *a, const blas_int *lda, const scs_float *beta, scs_float *c, const blas_int *ldc); +blas_int BLASC(herk)(const char *uplo, const char *trans, const blas_int *n, + const blas_int *k, const scs_float *alpha, + const scs_complex_float *a, const blas_int *lda, + const scs_float *beta, scs_complex_float *c, const blas_int *ldc); + void BLAS(scal)(const blas_int *n, const scs_float *sa, scs_float *sx, const blas_int *incx); +void BLASC(scal)(const blas_int *n, const scs_complex_float *sa, scs_complex_float *sx, + const blas_int *incx); #ifdef __cplusplus } @@ -48,6 +60,8 @@ void SCS(free_cone)(ScsCone *k) { scs_free(k->q); if (k->s) scs_free(k->s); + if (k->cs) + scs_free(k->cs); if (k->p) scs_free(k->p); scs_free(k); @@ -80,6 +94,13 @@ void SCS(deep_copy_cone)(ScsCone *dest, const ScsCone *src) { } else { dest->s = SCS_NULL; } + /* copy complex PSD cone */ + if (src->cssize > 0) { + dest->cs = (scs_int *)scs_calloc(src->cssize, sizeof(scs_int)); + memcpy(dest->cs, src->cs, src->cssize * sizeof(scs_int)); + } else { + dest->cs = SCS_NULL; + } /* copy power cone */ if (src->psize > 0) { dest->p = (scs_float *)scs_calloc(src->psize, sizeof(scs_float)); @@ -126,15 +147,19 @@ static inline scs_int get_sd_cone_size(scs_int s) { return (s * (s + 1)) / 2; } +static inline scs_int get_csd_cone_size(scs_int cs) { + return cs * cs; +} + /* * boundaries will contain array of indices of rows of A corresponding to * cone boundaries, boundaries[0] is starting index for cones of size strictly * larger than 1, boundaries malloc-ed here so should be freed. */ void set_cone_boundaries(const ScsCone *k, ScsConeWork *c) { - scs_int i, s_cone_sz, count = 0; + scs_int i, s_cone_sz, cs_cone_sz, count = 0; scs_int cone_boundaries_len = - 1 + k->qsize + k->ssize + k->ed + k->ep + k->psize; + 1 + k->qsize + k->ssize + k->cssize + k->ed + k->ep + k->psize; scs_int *b = (scs_int *)scs_calloc(cone_boundaries_len, sizeof(scs_int)); /* cones that can be scaled independently */ b[count] = k->z + k->l + k->bsize; @@ -147,7 +172,7 @@ void set_cone_boundaries(const ScsCone *k, ScsConeWork *c) { s_cone_sz = get_sd_cone_size(k->s[i]); b[count + i] = s_cone_sz; } - count += k->ssize; /* add ssize here not ssize * (ssize + 1) / 2 */ + count += k->ssize; /* size here not ssize * (ssize + 1) / 2 */ /* exp cones */ for (i = 0; i < k->ep + k->ed; ++i) { b[count + i] = 3; @@ -158,6 +183,11 @@ void set_cone_boundaries(const ScsCone *k, ScsConeWork *c) { b[count + i] = 3; } count += k->psize; + for (i = 0; i < k->cssize; ++i) { /*adding the complex psd cone here for backwards compatibility */ + cs_cone_sz = get_csd_cone_size(k->cs[i]); + b[count + i] = cs_cone_sz; + } + count += k->cssize; /* other cones */ c->cone_boundaries = b; c->cone_boundaries_len = cone_boundaries_len; @@ -175,6 +205,11 @@ static scs_int get_full_cone_dims(const ScsCone *k) { c += get_sd_cone_size(k->s[i]); } } + if (k->cssize) { + for (i = 0; i < k->cssize; ++i) { + c += get_csd_cone_size(k->cs[i]); + } + } if (k->ed) { c += 3 * k->ed; } @@ -238,6 +273,18 @@ scs_int SCS(validate_cones)(const ScsData *d, const ScsCone *k) { } } } + if (k->cssize && k->cs) { + if (k->cssize < 0) { + scs_printf("complex psd cone dimension error\n"); + return -1; + } + for (i = 0; i < k->cssize; ++i) { + if (k->cs[i] < 0) { + scs_printf("complex psd cone dimension error\n"); + return -1; + } + } + } if (k->ed && k->ed < 0) { scs_printf("ep cone dimension error\n"); return -1; @@ -266,15 +313,33 @@ void SCS(finish_cone)(ScsConeWork *c) { if (c->Xs) { scs_free(c->Xs); } + if (c->cXs) { + scs_free(c->cXs); + } if (c->Z) { scs_free(c->Z); } + if (c->cZ) { + scs_free(c->cZ); + } if (c->e) { scs_free(c->e); } + if (c->isuppz) { + scs_free(c->isuppz); + } if (c->work) { scs_free(c->work); } + if (c->iwork) { + scs_free(c->iwork); + } + if (c->cwork) { + scs_free(c->cwork); + } + if (c->rwork) { + scs_free(c->rwork); + } #endif if (c->cone_boundaries) { scs_free(c->cone_boundaries); @@ -289,7 +354,7 @@ void SCS(finish_cone)(ScsConeWork *c) { char *SCS(get_cone_header)(const ScsCone *k) { char *tmp = (char *)scs_malloc(512); /* sizeof(char) = 1 */ - scs_int i, soc_vars, sd_vars; + scs_int i, soc_vars, sd_vars, csd_vars; sprintf(tmp, "cones: "); if (k->z) { sprintf(tmp + strlen(tmp), "\t z: primal zero / dual free vars: %li\n", @@ -317,6 +382,14 @@ char *SCS(get_cone_header)(const ScsCone *k) { sprintf(tmp + strlen(tmp), "\t s: psd vars: %li, ssize: %li\n", (long)sd_vars, (long)k->ssize); } + csd_vars = 0; + if (k->cssize && k->cs) { + for (i = 0; i < k->cssize; i++) { + csd_vars += get_csd_cone_size(k->cs[i]); + } + sprintf(tmp + strlen(tmp), "\t cs: complex psd vars: %li, ssize: %li\n", + (long)csd_vars, (long)k->cssize); + } if (k->ep || k->ed) { sprintf(tmp + strlen(tmp), "\t e: exp vars: %li, dual exp vars: %li\n", (long)(3 * k->ep), (long)(3 * k->ed)); @@ -414,7 +487,7 @@ static scs_int proj_semi_definite_cone(scs_float *X, const scs_int n, #ifdef USE_LAPACK - /* copy lower triangular matrix into full matrix */ + /* copy upper triangular matrix into full matrix */ for (i = 0; i < n; ++i) { memcpy(&(Xs[i * (n + 1)]), &(X[i * n - ((i - 1) * i) / 2]), (n - i) * sizeof(scs_float)); @@ -466,7 +539,7 @@ static scs_int proj_semi_definite_cone(scs_float *X, const scs_int n, /* undo rescaling: scale diags by 1/sqrt(2) */ BLAS(scal)(&nb, &sqrt2_inv, Xs, &nb_plus_one); /* not n_squared */ - /* extract just lower triangular matrix */ + /* extract just upper triangular matrix */ for (i = 0; i < n; ++i) { memcpy(&(X[i * n - ((i - 1) * i) / 2]), &(Xs[i * (n + 1)]), (n - i) * sizeof(scs_float)); @@ -476,7 +549,171 @@ static scs_int proj_semi_definite_cone(scs_float *X, const scs_int n, #else scs_printf("FAILURE: solving SDP but no blas/lapack libraries were found!\n"); scs_printf("SCS will return nonsense!\n"); - SCS(scale_array)(X, NAN, n); + SCS(scale_array)(X, NAN, get_sd_cone_size(n)); + return -1; +#endif +} + +static scs_int set_up_csd_cone_work_space(ScsConeWork *c, const ScsCone *k) { + scs_int i; +#ifdef USE_LAPACK + blas_int n_max = 1; + blas_int neg_one = -1; + blas_int info = 0; + scs_float abstol = -1.0; + blas_int m = 0; + scs_complex_float lcwork = 0.0; + scs_float lrwork = 0.0; + blas_int liwork = 0; +#if VERBOSITY > 0 +#define _STR_EXPAND(tok) #tok +#define _STR(tok) _STR_EXPAND(tok) + scs_printf("BLAS(func) = '%s'\n", _STR(BLASC(func))); +#endif + /* eigenvector decomp workspace */ + for (i = 0; i < k->cssize; ++i) { + if (k->cs[i] > n_max) { + n_max = (blas_int)k->cs[i]; + } + } + c->cXs = (scs_complex_float *)scs_calloc(n_max * n_max, sizeof(scs_complex_float)); + c->cZ = (scs_complex_float *)scs_calloc(n_max * n_max, sizeof(scs_complex_float)); + c->e = (scs_float *)scs_calloc(n_max, sizeof(scs_float)); + c->isuppz = (blas_int *)scs_calloc(MAX(2, 2 * n_max), sizeof(blas_int)); + + /* workspace query */ + BLASC(heevr)("V", "A", "L", &n_max, c->cXs, &n_max, SCS_NULL, SCS_NULL, SCS_NULL, SCS_NULL, &abstol, &m, c->e, + c->cZ, &n_max, c->isuppz, &lcwork, &neg_one, &lrwork, &neg_one, &liwork, &neg_one, &info); + + if (info != 0) { + scs_printf("FATAL: heev workspace query failure, info = %li\n", (long)info); + return -1; + } + c->lcwork = (blas_int)(lcwork); + c->lrwork = (blas_int)(lrwork); + c->liwork = liwork; + c->cwork = (scs_complex_float *)scs_calloc(c->lcwork, sizeof(scs_complex_float)); + c->rwork = (scs_float *)scs_calloc(c->lrwork, sizeof(scs_float)); + c->iwork = (blas_int *)scs_calloc(c->liwork, sizeof(blas_int)); + + if (!c->cXs || !c->cZ || !c->e || !c->isuppz || !c->cwork || !c->rwork || !c->iwork) { + return -1; + } + return 0; +#else + for (i = 0; i < k->cssize; i++) { + if (k->cs[i] > 1) { + scs_printf( + "FATAL: Cannot solve SDPs without linked blas+lapack libraries\n"); + scs_printf( + "Install blas+lapack and re-compile SCS with blas+lapack library " + "locations\n"); + return -1; + } + } + return 0; +#endif +} + +/* size of X is get_csd_cone_size(n) */ +static scs_int proj_complex_semi_definite_cone(scs_float *X, const scs_int n, + ScsConeWork *c) { +/* project onto the positive semi-definite cone */ +#ifdef USE_LAPACK + scs_int i, first_idx; + blas_int nb = (blas_int)n; + blas_int ncols_z; + blas_int nb_plus_one = (blas_int)(n + 1); + blas_int one_int = 1; + scs_float zero = 0., one = 1.; + scs_complex_float csqrt2 = SQRTF(2.0); + scs_complex_float csqrt2_inv = 1.0 / csqrt2; + scs_complex_float *cXs = c->cXs; + scs_complex_float *cZ = c->cZ; + scs_float *e = c->e; + scs_float abstol = -1.0; + blas_int m = 0; + blas_int info = 0; + scs_complex_float csq_eig_pos; + +#endif + + if (n == 0) { + return 0; + } + if (n == 1) { + X[0] = MAX(X[0], 0.); + return 0; + } + +#ifdef USE_LAPACK + + /* copy upper triangular matrix into full matrix */ + for (i = 0; i < n - 1; ++i) { + cXs[i * (n + 1)] = X[i * (2 * n - i)]; + memcpy(&(cXs[i * (n + 1) + 1]), &(X[i * (2 * n - i) + 1]), + 2 * (n - i - 1) * sizeof(scs_float)); + } + cXs[n * n - 1] = X[n * n - 1]; + /* + rescale so projection works, and matrix norm preserved + see http://www.seas.ucla.edu/~vandenbe/publications/mlbook.pdf pg 3 + */ + /* scale diags by sqrt(2) */ + BLASC(scal)(&nb, &csqrt2, cXs, &nb_plus_one); /* not n_squared */ + + /* Solve eigenproblem, reuse workspaces */ + BLASC(heevr)("V", "A", "L", &nb, cXs, &nb, SCS_NULL, SCS_NULL, SCS_NULL, SCS_NULL, &abstol, &m, e, + cZ, &nb, c->isuppz, c->cwork, &c->lcwork, c->rwork, &c->lrwork, c->iwork, &c->liwork, &info); + if (info != 0) { + scs_printf("WARN: LAPACK heev error, info = %i\n", (int)info); + if (info < 0) { + return info; + } + } + + first_idx = -1; + /* e is eigvals in ascending order, find first entry > 0 */ + for (i = 0; i < n; ++i) { + if (e[i] > 0) { + first_idx = i; + break; + } + } + + if (first_idx == -1) { + /* there are no positive eigenvalues, set X to 0 and return */ + memset(X, 0, sizeof(scs_float) * get_csd_cone_size(n)); + return 0; + } + + /*LAPACK is col-major, so the columns of cZ' are the eigenvectors */ + /* scale cZ by sqrt(eig) */ + for (i = first_idx; i < n; ++i) { + csq_eig_pos = SQRTF(e[i]); + BLASC(scal)(&nb, &csq_eig_pos, &cZ[i * n], &one_int); + } + + /* Xs = Z Z' = V E V' */ + ncols_z = (blas_int)(n - first_idx); + BLASC(herk)("Lower", "NoTrans", &nb, &ncols_z, &one, &cZ[first_idx * n], &nb, &zero, cXs, &nb); + + /* undo rescaling: scale diags by 1/sqrt(2) */ + BLASC(scal)(&nb, &csqrt2_inv, cXs, &nb_plus_one); /* not n_squared */ + + /* extract just upper triangular matrix */ + for (i = 0; i < n - 1; ++i) { + X[i * (2 * n - i)] = cXs[i * (n + 1)]; + memcpy(&(X[i * (2 * n - i) + 1]), &(cXs[i * (n + 1) + 1]), + 2 * (n - i - 1) * sizeof(scs_float)); + } + X[n * n - 1] = cXs[n * n - 1]; + return 0; + +#else + scs_printf("FAILURE: solving SDP but no blas/lapack libraries were found!\n"); + scs_printf("SCS will return nonsense!\n"); + SCS(scale_array)(X, NAN, get_csd_cone_size(n)); return -1; #endif } @@ -720,6 +957,17 @@ static scs_int proj_cone(scs_float *x, const ScsCone *k, ScsConeWork *c, } } + if (k->cssize && k->cs) { /* doesn't use r_y */ + /* project onto complex PSD cones */ + for (i = 0; i < k->cssize; ++i) { + status = proj_complex_semi_definite_cone(&(x[count]), k->cs[i], c); + if (status < 0) { + return status; + } + count += get_csd_cone_size(k->cs[i]); + } + } + if (k->ep || k->ed) { /* doesn't use r_y */ #ifdef _OPENMP #pragma omp parallel for @@ -775,6 +1023,12 @@ ScsConeWork *SCS(init_cone)(ScsCone *k, scs_int m) { return SCS_NULL; } } + if (k->cssize && k->cs) { + if (set_up_csd_cone_work_space(c, k) < 0) { + SCS(finish_cone)(c); + return SCS_NULL; + } + } return c; }