26static int dsdpuselapack=1;
28#define __FUNCT__ "DSDPUseLAPACKForSchur"
29int DSDPUseLAPACKForSchur(
DSDP dsdp,
int flag){
32 DSDPFunctionReturn(0);
36#define __FUNCT__ "Taddline"
37static int Taddline(
void* M,
int row,
double dd,
double v[],
int m){
39 MCholSolverALL *AMA = (MCholSolverALL *)M;
41 info=MatAddColumn4(AMA->M,dd,v,row); DSDPCHKERR(info);
42 DSDPFunctionReturn(0);
46#define __FUNCT__ "Tadddiagonal"
47static int Tadddiagonal(
void* M,
int row,
double v){
49 MCholSolverALL *AMA = (MCholSolverALL *)M;
51 info=MatAddDiagonalElement(AMA->M,row,v); DSDPCHKERR(info);
52 DSDPFunctionReturn(0);
56#define __FUNCT__ "Tassemble"
57static int Tassemble(
void*M){
59 DSDPFunctionReturn(0);
64#define __FUNCT__ "DSDPCheckForSparsity"
65static int DSDPCheckForSparsity(
DSDP dsdp,
int m,
int *rnnz,
int tnnz[],
int *totalnnz){
66 int info,i,j,tottalnnz=0;
68 memset(rnnz,0,(m+1)*
sizeof(
int));
71 for (j=i+1;j<m;j++){
if (tnnz[j]>0) rnnz[i+1]++;}
74 for (i=0;i<m;i++){tottalnnz+=rnnz[i+1];}
76 DSDPFunctionReturn(0);
80#define __FUNCT__ "DSDPCreateM"
81static int DSDPCreateM(MCholSolverALL *ABA, chfac **M,
int rrnnz[],
int tnnz[],
int totalnnz){
87 int n=ABA->m,m=ABA->m;
92 DSDPCALLOC2(&snnz,
int,(totalnnz+1),&info); DSDPCHKERR(info);
93 DSDPCALLOC2(&rnnz,
int,(m+1),&info); DSDPCHKERR(info);
94 for (i=0;i<=m;i++){ rnnz[i]=rrnnz[i];}
98 for (j=i+1;j<m;j++){
if (tnnz[j]>0) {snnz[tt]=j; tt++;} }
101 printf(
"Trying Sparse M: Total nonzeros: %d of %d \n",totalnnz,m*(m-1)/2 );
103 SymbProc(rnnz+1,snnz,n,&sfptr);
116 for (i=m-1;i>=0;i--){
117 for (j=rnnz[i+1]-1;j>=rnnz[i];j--){
121 for (k=j;k<rnnz[col]-1;k++){ snnz[k]=snnz[k+1];}
122 for (k=i;k<col;k++){ rnnz[k+1]--;}
128 DSDPFunctionReturn(0);
133#define __FUNCT__ "DSDPLinearSolverPrepare"
134static int DSDPLinearSolverPrepare(
void* ctx,
int*flag){
138 MCholSolverALL *AMA = (MCholSolverALL *)ctx;
146 Cfact=(cfc_sta)ChlFact(sf,sf->iw,sf->rw,TRUE);
147 if (CfcOk!=Cfact ){ *flag=1; }
148 DSDPFunctionReturn(0);
153#define __FUNCT__ "DSDPLinearSolve"
154static int DSDPLinearSolve2(
void* ctx,
double dd[],
double xx[],
int n){
158 MCholSolverALL *AMA = (MCholSolverALL *)ctx;
161 info=DSDPVecGetArray(AMA->D1, &bb);DSDPCHKERR(info);
162 for (i=0;i<n;i++){ bb[i]=dd[i];}
163 ChlSolve(AMA->M, bb, xx);
164 info=DSDPVecRestoreArray(AMA->D1, &bb);
165 DSDPFunctionReturn(0);
170#define __FUNCT__ "DSDPGramMatRowNonzeros"
171static int DSDPGramMatRowNonzeros(
void*M,
int row,
double cols[],
int *ncols,
int nrows){
172 MCholSolverALL *AMA = (MCholSolverALL *)M;
177 for (i=row;i<nrows;i++){
181 *ncols = AMA->rnnz[row+1] - AMA->rnnz[row]+1;
183 for (i=AMA->rnnz[row]; i<AMA->rnnz[row+1]; i++){
184 cols[AMA->colnnz[i]]=1.0;
187 DSDPFunctionReturn(0);
191#define __FUNCT__ "Tzero"
192static int Tzero(
void*M){
194 MCholSolverALL *AMA = (MCholSolverALL *)M;
196 info=MatZeroEntries4(AMA->M); DSDPCHKERR(info);
197 DSDPFunctionReturn(0);
201#define __FUNCT__ "Tdestroy"
202static int Tdestroy(
void*M){
203 MCholSolverALL *AMA = (MCholSolverALL *)M;
207 info = DSDPVecDestroy(&AMA->D1); DSDPCHKERR(info);
208 if (AMA->isdense==0 && AMA->rnnz ){
209 DSDPFREE(&AMA->rnnz,&info);DSDPCHKERR(info);
210 DSDPFREE(&AMA->colnnz,&info);DSDPCHKERR(info);
212 DSDPFREE(&AMA,&info);DSDPCHKERR(info);
213 DSDPFunctionReturn(0);
216static int Tsetup(
void*M,
int m){
218 DSDPFunctionReturn(0);
221static int TTTMatMult(
void*M,
double x[],
double y[],
int n){
222 MCholSolverALL *AMA = (MCholSolverALL *)M;
225 info=MatMult4(AMA->M,x,y,n); DSDPCHKERR(info);
226 DSDPFunctionReturn(0);
229static int TTTMatShiftDiagonal(
void *M,
double d){
230 MCholSolverALL *AMA = (MCholSolverALL *)M;
233 info=Mat4DiagonalShift(AMA->M,d); DSDPCHKERR(info);
234 DSDPFunctionReturn(0);
237static int TTTMatAddDiagonal(
void *M,
double diag[],
int m){
238 MCholSolverALL *AMA = (MCholSolverALL *)M;
241 info=Mat4AddDiagonal(AMA->M,diag,m); DSDPCHKERR(info);
242 DSDPFunctionReturn(0);
245static int TTTMatView(
void *M){
246 MCholSolverALL *AMA = (MCholSolverALL *)M;
249 info=Mat4View(AMA->M); DSDPCHKERR(info);
250 DSDPFunctionReturn(0);
263static const char* tmatname=
"SPARSE PSD";
265static int TMatOpsInit(
struct DSDPSchurMat_Ops *mops){
269 mops->matrownonzeros=DSDPGramMatRowNonzeros;
270 mops->mataddrow=Taddline;
271 mops->matadddiagonal=TTTMatAddDiagonal;
272 mops->mataddelement=Tadddiagonal;
273 mops->matshiftdiagonal=TTTMatShiftDiagonal;
274 mops->matassemble=Tassemble;
275 mops->matscaledmultiply=TTTMatMult;
276 mops->matfactor=DSDPLinearSolverPrepare;
277 mops->matsolve=DSDPLinearSolve2;
278 mops->matdestroy=Tdestroy;
280 mops->matsetup=Tsetup;
281 mops->matview=TTTMatView;
283 mops->matname=tmatname;
284 DSDPFunctionReturn(0);
287static struct DSDPSchurMat_Ops dsdpmatops;
289int DSDPGetDiagSchurMat(
int,
struct DSDPSchurMat_Ops **,
void **);
290int DSDPGetLAPACKSUSchurOps(
int,
struct DSDPSchurMat_Ops**,
void**);
291int DSDPGetLAPACKPUSchurOps(
int,
struct DSDPSchurMat_Ops**,
void**);
293static int DSDPCreateM(MCholSolverALL*,chfac**,
int[],
int[],
int);
294static int DSDPCreateSchurMatrix(
void*,
int);
297#define __FUNCT__ "DSDPCreateSchurMatrix"
298static int DSDPCreateSchurMatrix(
void *ctx,
int m){
301 int *rnnz,*tnnz,totalnnz;
307 struct DSDPSchurMat_Ops *tmatops;
311 info=DSDPGetDiagSchurMat(m,&tmatops,&tdata); DSDPCHKERR(info);
316 DSDPCALLOC2(&rnnz,
int,(m+1),&info); DSDPCHKERR(info);
317 DSDPCALLOC2(&tnnz,
int,(m+1),&info); DSDPCHKERR(info);
319 info=DSDPCheckForSparsity(dsdp,m,rnnz,tnnz,&totalnnz); DSDPCHKERR(info);
321 if (totalnnz*2+m > m*m*0.1 && dsdpuselapack) {
323 info=DSDPGetLAPACKSUSchurOps(m,&tmatops,&tdata);
324 if (info) {gotit=0;printf(
"Try packed format\n"); }
325 DSDPLogInfo(0,8,
"Creating Dense Full LAPACK Schur Matrix\n");
329 if ( 0 && totalnnz*2+m > m*m*0.1 && dsdpuselapack) {
331 info=DSDPGetLAPACKPUSchurOps(m,&tmatops,&tdata); DSDPCHKERR(info);
332 DSDPLogInfo(0,8,
"Creating Dense Packed LAPACK Schur Matrix\n");
339 DSDPCALLOC1(&AMA,MCholSolverALL,&info);DSDPCHKERR(info);
340 AMA->dsdp=dsdp; AMA->m=m;
341 info=DSDPVecCreateSeq(m,&AMA->D1); DSDPCHKERR(info);
342 if (totalnnz*2+m > m*m * 0.11 ){
343 info=MchlSetup2(m,&sf); DSDPCHKERR(info);
344 AMA->M=sf; AMA->isdense=1;
345 AMA->rnnz=0; AMA->colnnz=0;
346 DSDPLogInfo(0,8,
"Creating Dense Full non LAPACK Schur Matrix\n");
348 info=DSDPCreateM(AMA,&sf,rnnz,tnnz,totalnnz); DSDPCHKERR(info);
349 DSDPLogInfo(0,8,
"Creating Sparse Schur Matrix\n");
352 info=TMatOpsInit(&dsdpmatops); DSDPCHKERR(info);
355 DSDPFREE(&tnnz,&info);DSDPCHKERR(info);
356 DSDPFREE(&rnnz,&info);DSDPCHKERR(info);
357 DSDPFunctionReturn(0);
360static struct DSDPSchurMat_Ops dsdpmatops000;
363#define __FUNCT__ "DSDPSetDefaultSchurMatrixStructure"
364int DSDPSetDefaultSchurMatrixStructure(
DSDP dsdp){
368 dsdpmatops000.matsetup=DSDPCreateSchurMatrix;
370 DSDPFunctionReturn(0);
int DSDPSetSchurMatOps(DSDP, struct DSDPSchurMat_Ops *, void *)
Set the Schur complement matrix.
The API to DSDP for those applications using DSDP as a subroutine library.
struct DSDP_C * DSDP
An implementation of the dual-scaling algorithm for semidefinite programming.
int DSDPSparsityInSchurMat(DSDP dsdp, int row, int rnnz[], int mm)
Identify nonzero elements in a row of the Schur complement.
int DSDPSchurMatOpsInitialize(struct DSDPSchurMat_Ops *dops)
Initialize function pointers to 0.
Function pointers that a Schur complement matrix (dense, sparse, parallel dense) must provide.
Error handling, printing, and profiling.
Vector operations used by the solver.
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Internal structures for the DSDP solver.