/
blksdp.h
210 lines (192 loc) · 8.43 KB
/
blksdp.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
/* ************************************************************
HEADER blksdp.h
For use with mex-files in self-dual-minimization package.
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
% Dept. Econometrics & O.R., Tilburg University, the Netherlands.
% Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
% Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
% CRL, McMaster University, Canada.
% Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
************************************************************ */
#if !defined(BLKSDP)
#define BLKSDP
#include "mex.h"
#ifdef OCTAVE
#include "f77blas.h" /* defines "blasint" data type */
#define FORT(x) BLASFUNC(x)
#else /* Matlab */
#include "blas.h"
typedef ptrdiff_t blasint;
/**
* For Matlab R2019a (probably before) and newer, when including
* "blas.h" the respective BLAS identifiers are already defined,
* e.g. for "dcopy":
*
* #define dcopy FORTRAN_WRAPPER(dcopy)
*
* thus calling FORT(dcopy) == FORTRAN_WRAPPER(dcopy) inside SeDuMi
* would result in "dcopy__" and already resulted in some bug reports.
*
* Compiling with -DFWRAPPER restores the previous behavior.
*/
#ifdef FWRAPPER
#define FORT(x) FORTRAN_WRAPPER(x)
#else
#define FORT(x) x
#endif
#endif
/* ------------------------------------------------------------
Type definitions:
------------------------------------------------------------ */
typedef struct{
double *pr;
mwIndex *jc, *ir;
} jcir;
typedef struct{
double *pr;
const mwIndex *jc;
const mwIndex *ir;
} constjcir;
typedef struct{
mwSize frN,lpN,lorN,rconeN,sdpN, rsdpN;
mwSize qMaxn,rMaxn,hMaxn, rLen,hLen, qDim,rDim,hDim;
const double *lorNL,*rconeNL,*sdpNL;
} coneK;
/* ------------------------------------------------------------
Macros:
------------------------------------------------------------ */
#if !defined(SQR)
#define SQR(x) ((x)*(x))
#endif
#if !defined(MAX)
#define MAX(A, B) ((A) > (B) ? (A) : (B))
#endif
#if !defined(MIN)
#define MIN(A, B) ((A) < (B) ? (A) : (B))
#endif
#if !defined(SIGN)
#define SIGN(A) (2 * ((A) >= 0) - 1)
#endif
#ifndef M_SQRT2
#define M_SQRT2 1.41421356237309504880 /* sqrt(2) */
#endif
#ifndef M_SQRT1_2
#define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */
#endif
/* ************************************************************
INT COMPARE: for searching an mwIndex array
NOTE: qsort sorts in ascending (0,1,2,..) order, if the compare
function returns < 0 iff a<b, 0 iff a==b, > 0 iff a > b.
************************************************************ */
#if !defined(_COMPFUN_)
#define _COMPFUN_
typedef int (*COMPFUN)(const void *pa,const void *pb);
#endif
#define ibsearch(key,vec,n) bsearch((void *)(key), (void *)(vec), (n), sizeof(mwSize), (COMPFUN) icmp)
#define iqsort(vec,n) qsort((void *)(vec), (n), sizeof(mwSize), (COMPFUN) icmp)
/* --------------------------------------
KEY COMPARE: FOR SORTING AN (INT or FLOAT) ARRAY WITH INT-KEYS.
-------------------------------------- */
typedef struct{
mwIndex i,k;
} keyint;
#if !defined(_KEYDOUBLE_)
#define _KEYDOUBLE
typedef struct{
double r;
mwIndex k;
} keydouble;
#endif
#define kiqsort(vec,n) qsort((void *)(vec), (n), sizeof(keyint), (COMPFUN) kicmp);
#define kdsortdec(vec,n) qsort((void *)(vec), (n), sizeof(keydouble), (COMPFUN) kdcmpdec);
/* ------------------------------------------------------------
Prototypes:
------------------------------------------------------------ */
char icmp(const mwIndex *a, const mwIndex *b);
bool intbsearch(mwIndex *pi, const mwIndex *x, const mwIndex n, const mwIndex key);
char intmbsearch(mwIndex *z, bool *found, const mwIndex *x, const mwIndex xnnz,
const mwIndex *y, const mwIndex ynnz, mwIndex *iwork, const mwIndex iwsize);
char kicmp(const keyint *a, const keyint *b);
char kdcmpdec(const keydouble *a, const keydouble *b);
double realssqr(const double *x, const mwSize n);
double realdot(const double *x, const double *y, const mwSize n);
double selrealdot(const double *x, const double *y,
const mwIndex *sel, const mwSize nnz);
double realdotrow(const double *x, const double *y, const mwSize n);
void fromto(mwIndex *x, mwSize i, const mwSize n);
double triudotprod(const double *x, const double *y, const mwSize n);
double striudotprod(const double *x, const double *y, const mwSize n);
void tril2sym(double *r, const mwSize n);
void tril2herm(double *r, double *rpi, const mwSize n);
void triu2sym(double *r, const mwSize n);
void triu2herm(double *r, double *rpi, const mwSize n);
void scalarmul(double *r, const double alpha,const double *x,const mwSize n);
void addscalarmul(double *r, const double alpha,const double *x,const mwSize n);
void subscalarmul(double *x, const double alpha, const double *y, const mwSize n);
void realHadamard(double * r, const double *x, const double *y, const mwSize n);
void minusHadamard(double * r, const double *x, const double *y, const mwSize n);
void realHadarow(double * r, const double *x, const double *y, const mwSize n);
void realHadadiv(double * r, const double *x, const double *y, const mwSize n);
void fzeros(double *z,const mwSize n);
void conepars(const mxArray *mxK, coneK *pK);
void someStats(mwSize *pxmax, mwIndex *pxsum, mwIndex *pxssqr,
const double *x, const mwSize n);
mwIndex spsqrscale(double *z, mwIndex *blks, const mwIndex *zjc, const mwIndex *zir,
const mwIndex *znnz, const double *d,
const mwIndex *xir, const double *xpr, mwIndex xjc0, const mwIndex xjc1,
const mwIndex *blkstart, const mwIndex *xblk, const mwIndex *psdNL,
const mwIndex rpsdN, double *fwork, mwIndex *iwork);
#ifdef OLDSEDUMI
double qscale(double *z,const double *x,const double *y,
const double rdetx,const mwIndex n);
void qlmul(double *z,const double *x,const double *y,
const double rdetx,const mwIndex n);
void qldiv(double *z,const double *x,const double *y,
const double rdetx,const mwIndex n);
void vec2blks(mwIndex *blklocs, const mwIndex *blkstart, const mwIndex *yir,
const mwIndex ystart, const mwIndex ynnz, const mwIndex nblk);
void vec2selblks(mwIndex *blklocs, const mwIndex *blkstart, const mwIndex *yir,
const mwIndex ystart, const mwIndex ynnz,
const mwIndex *blkir, const mwIndex blknnz);
mwIndex lqdsqrx(double *z,
const mwIndex *xir, const double *xpr, const mwIndex xjc0,
const mwIndex xjcq, const mwIndex xjcs, const mwIndex *qir,
const mwIndex *blkstart,
const double *dsqr, const double *detd);
mwIndex blkpsdscale(double *z, const mwIndex *zir, const mwIndex zjc1,
const double *u, const mwIndex *invperm, const double *x,
const mwIndex *xblk, const mwIndex blkjc0, const mwIndex blkjc1,
const mwIndex *blkstart, const mwIndex *psdNL, const mwIndex *cumpsdNL,
const mwIndex rpsdN, double *fwork);
#endif
void uperm(double *y, const double *u, const mwIndex *perm, const mwIndex n);
/* ------------------------------------------------------------
For auxfwdpr1:
------------------------------------------------------------ */
void fwipr1(double *y, const double *p, const double *beta,
const mwSize m, const mwSize n);
void fwipr1o(double *y, const mwIndex *perm, const double *p, const double *beta,
const mwSize m, const mwSize n);
#endif