My Project
OSCsdpSolver.cpp
Go to the documentation of this file.
1/* $Id: OSCsdpSolver.cpp 4787 2014-01-22 22:26:56Z Gassmann $ */
20#include "OSCsdpSolver.h"
21#include "OSFileUtil.h"
22#include "OSInstance.h"
23#include "OSOption.h"
24#include "OSGeneral.h"
25#include "OSMatrix.h"
26#include "OSrLReader.h"
27#include "OSOutput.h"
28#include "OSParameters.h"
29#include "OSMathUtil.h"
30
31#include "CoinTime.hpp"
32
33#include <map>
34
35#include <iostream>
36#ifdef HAVE_CTIME
37# include <ctime>
38#else
39# ifdef HAVE_TIME_H
40# include <time.h>
41# else
42# error "don't have header file for time"
43# endif
44#endif
45
46using std::cout;
47using std::endl;
48using std::ostringstream;
49
50
52{
53#ifndef NDEBUG
55 "inside CsdpSolver constructor\n");
56#endif
57 osrlwriter = new OSrLWriter();
58 osresult = new OSResult();
59 m_osilreader = NULL;
60 m_osolreader = NULL;
61 csdpErrorMsg = new std::string("");
62}
63
65{
66#ifndef NDEBUG
67 osoutput->OSPrint(ENUM_OUTPUT_AREA_OSSolverInterfaces, ENUM_OUTPUT_LEVEL_debug, "inside CsdpSolver destructor\n");
68#endif
69 if(m_osilreader != NULL) delete m_osilreader;
70 m_osilreader = NULL;
71 if(m_osolreader != NULL) delete m_osolreader;
72 m_osolreader = NULL;
73 delete osresult;
74 osresult = NULL;
75 delete osrlwriter;
76 osrlwriter = NULL;
77 //delete osinstance;
78 //osinstance = NULL;
79 delete csdpErrorMsg;
80#ifndef NDEBUG
81 osoutput->OSPrint(ENUM_OUTPUT_AREA_OSSolverInterfaces, ENUM_OUTPUT_LEVEL_trace, "Leaving CsdpSolver destructor\n");
82#endif
83}
84
86{
87 std::ostringstream outStr;
88 ScalarExpressionTree* tempTree;
89 OSnLNode *tr;
90 OSnLMNode *mt;
91 OSnLMNode *mr;
92 OSnLMNode *mv;
93 OSMatrix* tempMtx;
94
95 int* blockOffset = NULL;
96 int* blockSize = NULL;
97 int* mtxRef = NULL;
98 bool* isdiag = NULL;
99 ExpandedMatrixBlocks** mtxBlocks = NULL;
100 try
101 {
102 if(osil.length() == 0 && osinstance == NULL) throw ErrorClass("there is no instance");
103 clock_t start, finish;
104 double duration;
105 start = clock();
106 if(osinstance == NULL)
107 {
108 m_osilreader = new OSiLReader();
110 }
111 finish = clock();
112 duration = (double) (finish - start) / CLOCKS_PER_SEC;
113
114 /* Process the osinstance into the --- rather tricky --- CSDP data structures
115 * and verify that the solver is appropriate - CSDP requires a very special type of problem
116 */
117
118 // Check general problem characteristics
120 throw ErrorClass("There must be one matrixVar object");
122 throw ErrorClass("There must be one nonlinear expression for each constraint and objective");
124 throw ErrorClass("Additional linear constraint coefficients are not supported");
126 throw ErrorClass("Additional quadratic terms are not supported");
127
128 char* cType = osinstance->getConstraintTypes();
129 for (int i=0; i < osinstance->getConstraintNumber(); i++)
130 if (cType[i] != 'E') throw ErrorClass("Only equality constraints are supported");
131
132 std::string* oType = osinstance->getObjectiveMaxOrMins();
133 if (oType[0] != "max") throw ErrorClass("The problem must be of \"max\" type");
134
135 //Check the form of the objective
137 if (tempTree == NULL) throw ErrorClass("Expecting matrixTrace in objective row");
138 tr = tempTree->m_treeRoot;
139 if (tr->inodeInt != OS_MATRIX_TRACE)
140 throw ErrorClass("Expecting matrixTrace in objective row");
141 mt = tr->m_mMatrixChildren[0];
142 if (mt->inodeInt != OS_MATRIX_TIMES)
143 throw ErrorClass("Unsupported expression in objective row");
144 mr = mt->m_mMatrixChildren[0];
145 mv = mt->m_mMatrixChildren[1];
147 throw ErrorClass("Unsupported expression in objective row");
148
149 mtxRef = new int[osinstance->getConstraintNumber()+1];
150
151 // Analyze A0 matrix: Verify existence, block-diagonal structure, get block dimensions, etc.
152 mtxRef[0] = ((OSnLMNodeMatrixReference*)mr)->idx;
153 if (osinstance->instanceData->matrices == NULL)
154 throw ErrorClass("<matrices> section was never defined");
155 if (mtxRef[0] < 0 || mtxRef[0] >= osinstance->getMatrixNumber())
156 throw ErrorClass("Illegal matrix reference");
157 tempMtx = osinstance->instanceData->matrices->matrix[mtxRef[0]];
158 if (tempMtx == NULL) throw ErrorClass("A0 matrix was never defined");
159 if (tempMtx->numberOfRows != tempMtx->numberOfColumns)
160 throw ErrorClass("A0 matrix must be square and symmetric");
161 if (tempMtx->getMatrixType() != ENUM_MATRIX_TYPE_constant)
162 throw ErrorClass("A0 matrix must be of type \"constant\"");
163
164 int* rowOffsets = tempMtx->getRowPartition();
165 int nRowBlocks = tempMtx->getRowPartitionSize();
166 int* columnOffsets = tempMtx->getColumnPartition();
167 int nColumnBlocks = tempMtx->getColumnPartitionSize();
168
169 int* tempRowOffsets;
170 int tempNRowBlocks;
171 int* tempColumnOffsets;
172 int tempNColumnBlocks;
173
174 int i0, itemp, imerge;
175
176 //do the same for all constraints
177 for (int i=0; i < osinstance->getConstraintNumber(); i++)
178 {
180 if (tempTree == NULL) throw ErrorClass("Expecting matrixTrace in constraint row");
181 tr = tempTree->m_treeRoot;
182 if (tr->inodeInt != OS_MATRIX_TRACE)
183 throw ErrorClass("Expecting matrixTrace in constraint row");
184 mt = tr->m_mMatrixChildren[0];
185 if (mt->inodeInt != OS_MATRIX_TIMES)
186 throw ErrorClass("Unsupported expression in constraint row");
187 mr = mt->m_mMatrixChildren[0];
188 mv = mt->m_mMatrixChildren[1];
190 throw ErrorClass("Unsupported expression in constraint row");
191
192 // Analyze Ai matrix: Verify existence, block-diagonal structure, get block dimensions, etc.
193 mtxRef[i+1] = ((OSnLMNodeMatrixReference*)mr)->idx;
194 if (mtxRef[i+1] < 0 || mtxRef[i+1] >= osinstance->getMatrixNumber())
195 throw ErrorClass("Illegal matrix reference");
196 tempMtx = osinstance->instanceData->matrices->matrix[mtxRef[i+1]];
197 if (tempMtx == NULL) throw ErrorClass("Matrix in constraint was never defined");
198 if (tempMtx->numberOfRows != tempMtx->numberOfColumns)
199 throw ErrorClass("Constraint matrix must be square and symmetric");
200 if (tempMtx->getMatrixType() != ENUM_MATRIX_TYPE_constant)
201 throw ErrorClass("Constraint matrix must be of type \"constant\"");
202
203 tempRowOffsets = tempMtx->getRowPartition();
204 tempNRowBlocks = tempMtx->getRowPartitionSize();
205 tempColumnOffsets = tempMtx->getColumnPartition();
206 tempNColumnBlocks = tempMtx->getColumnPartitionSize();
207
208 // merge row partitions
209 i0 = 0;
210 itemp = 0;
211 imerge = 0;
212 for (;;)
213 {
214 if (rowOffsets[i0] == tempRowOffsets[itemp])
215 {
216 if (imerge != i0) rowOffsets[imerge] = rowOffsets[i0];
217 i0++;
218 itemp++;
219 imerge++;
220 }
221 else
222 {
223 if (rowOffsets[i0] < tempRowOffsets[itemp])
224 i0++;
225 else
226 itemp++;
227 }
228 if (i0 >= nRowBlocks || itemp >= tempNRowBlocks)
229 break;
230 }
231 nRowBlocks = imerge;
232
233 // merge column partititons
234 i0 = 0;
235 itemp = 0;
236 imerge = 0;
237 for (;;)
238 {
239 if (columnOffsets[i0] == tempColumnOffsets[itemp])
240 {
241 if (imerge != i0) columnOffsets[imerge] = columnOffsets[i0];
242 i0++;
243 itemp++;
244 imerge++;
245 }
246 else
247 {
248 if (columnOffsets[i0] < tempColumnOffsets[itemp])
249 i0++;
250 else
251 itemp++;
252 }
253 if (i0 >= nColumnBlocks || itemp >= tempNColumnBlocks)
254 break;
255 }
256 nColumnBlocks = imerge;
257 }
258
259 blockOffset = new int[nRowBlocks]; // step this through: nRowBlocks=?
260
261 // make sure the row and column blocks are synchronized and compute block sizes
262 int nBlocks;
263 int jrow = 0;
264 int jcol = 0;
265 nBlocks = 0;
266 for (;;)
267 {
268 if (rowOffsets[jrow] == columnOffsets[jcol])
269 {
270 blockOffset[nBlocks] = rowOffsets[jrow];
271 jrow++;
272 jcol++;
273 nBlocks++;
274 }
275 else
276 {
277 if (rowOffsets[jrow] < columnOffsets[jcol])
278 jrow++;
279 else
280 jcol++;
281 }
282 if (jrow >= nRowBlocks || jcol >= nColumnBlocks)
283 break;
284 }
285
286 // Note: nBlocks is one larger than the number of blocks.
287 // blockSize is 1-based, due to issues of Fortran compatibility,
288 // so we will burn off blockSize[0] anyway...
289 blockSize = new int[nBlocks];
290 for (int i=1; i < nBlocks; i++)
291 {
292#ifndef NOSHORTS
293 if (blockOffset[i] - blockOffset[i-1] >= USHRT_MAX)
294 throw ErrorClass("This problem is too large to be solved by this version of the code!\n"
295 + "Recompile without -DUSERSHORTINDS to fix the problem.\n");
296#endif
297 blockSize[i] = blockOffset[i] - blockOffset[i-1];
298 }
299
300 nC_rows = blockOffset[nBlocks-1];
301 nC_blks = nBlocks-1;
303
304#ifndef NOSHORTS
306 if (ncon >= USHRT_MAX)
307 throw ErrorClass("This problem is too large to be solved by this version of the code!\n"
308 + "Recompile without -DUSERSHORTINDS to fix the problem.\n");
309
310 if (nBlocks >= USHRT_MAX)
311 throw ErrorClass("This problem is too large to be solved by this version of the code!\n"
312 + "Recompile without -DUSERSHORTINDS to fix the problem.\n");
313#endif
314
315#ifndef BIT64
316 /*
317 * If operating in 32 bit mode, make sure that the dimension mDIM isn't
318 * too big for 32 bits. If we don't do this check, then integer overflow
319 * won't be detected, and we'll allocate a bogus amount of storage.
320 */
321 if (ncon > 23169)
322 throw ErrorClass("This problem is too large to be solved in 32 bit mode!\n");
323#endif
324
325 // set up the right hand side values. Note: 1-based for Fortran interface.
327
328 // Set up storage and retrieve pointers.
329 mtxBlocks = new ExpandedMatrixBlocks*[ncon+1];
330 GeneralSparseMatrix* tmpBlock;
331
332 // At this point we know the dimensions of all blocks.
333 // Keep track of diagonal blocks. Note: isdiag is 1-indexed
334 isdiag = new bool[nBlocks];
335 for (int i=1; i<nBlocks; i++)
336 isdiag[i] = true;
337
338 for (int j=0; j < ncon+1; j++)
339 {
340 mtxBlocks[j] = osinstance->instanceData->matrices->matrix[mtxRef[j]]
341 ->getBlocks(blockOffset,nBlocks,blockOffset,nBlocks,false,true);
342
343 if (!mtxBlocks[j]->isBlockDiagonal())
344 throw ErrorClass("Constraint matrix must be block-diagonal");
345
346 for (int i=0; i<nBlocks-1; i++)
347 {
348 tmpBlock = mtxBlocks[j]->getBlock(i,i);
349 if (tmpBlock != NULL && !(tmpBlock->isDiagonal()))
350 isdiag[i+1] = false;
351 }
352 }
353
355 C_matrix.nblocks=nBlocks-1;
356 C_matrix.blocks = new blockrec[nBlocks];
357
359 for (int blk=1; blk < nBlocks; blk++)
360 {
361 tmpBlock = mtxBlocks[0]->getBlock(blk-1,blk-1);
362 int blksz = blockSize[blk];
363 if (isdiag[blk] == 1)
364 {
365 // diagonal block
366 C_matrix.blocks[blk].blocksize = blksz;
367 C_matrix.blocks[blk].blockcategory = DIAG;
368 C_matrix.blocks[blk].data.vec = new double[blksz+1];
369
370 for (int i=1; i<=blksz; i++)
371 C_matrix.blocks[blk].data.vec[i] = 0.0;
372
373 if (tmpBlock != NULL)
374 {
375 for (int i=0; i < tmpBlock->valueSize; i++)
376 C_matrix.blocks[blk].data.vec[tmpBlock->index[i]+1]
377 = ((ConstantMatrixValues*)tmpBlock->value)->el[i];
378 }
379 }
380 else
381 {
382 // There are off-diagonals (i.e., "matrix block")
383 C_matrix.blocks[blk].blocksize = blksz;
384 C_matrix.blocks[blk].blockcategory = MATRIX;
385 C_matrix.blocks[blk].data.mat = new double[blksz*blksz];
386
387 for (int i=1; i<=blksz; i++)
388 for (int j=1; j<=blksz; j++)
389 C_matrix.blocks[blk].data.mat[ijtok(i,j,blksz)] = 0.0;
390
391 if (tmpBlock != NULL)
392 {
393 for (int i=1; i < tmpBlock->startSize; i++)
394 for (int j=tmpBlock->start[i-1]; j<tmpBlock->start[i]; j++)
395 {
396 C_matrix.blocks[blk].data.mat[ijtok(i,tmpBlock->index[j]+1,blksz)]
397 = ((ConstantMatrixValues*)tmpBlock->value)->el[j];
398 C_matrix.blocks[blk].data.mat[ijtok(tmpBlock->index[j]+1,i,blksz)]
399 = ((ConstantMatrixValues*)tmpBlock->value)->el[j];
400 }
401 }
402 }
403 }
404
406 mconstraints = new constraintmatrix[ncon+1];
407
409 for (int i=1; i<=ncon; i++)
410 {
411 mconstraints[i].blocks = NULL;
412 }
413
414 struct sparseblock *p;
415 struct sparseblock *q;
416 struct sparseblock *prev;
417
422 for (int i=1; i<=ncon; i++)
423 {
424 prev = NULL;
425 for (int blk=1; blk < nBlocks; blk++)
426 {
427 tmpBlock = mtxBlocks[i]->getBlock(blk-1,blk-1);
428 if (tmpBlock != NULL && tmpBlock->valueSize > 0)
429 {
434 p = new sparseblock();
435 p->numentries = tmpBlock->valueSize;
436 p->entries = new double[p->numentries+1];
437#ifdef NOSHORTS
438 p->iindices = new int[p->numentries+1];
439 p->jindices = new int[p->numentries+1];
440#else
441 p->iindices = new unsigned short[p->numentries+1];
442 p->jindices = new unsigned short[p->numentries+1];
443#endif
444 p->blocknum = blk;
445 p->blocksize = blockSize[blk];
446 p->constraintnum = i;
447 p->next = NULL;
448 p->nextbyblock = NULL;
449 if (((p->numentries) > 0.25*(p->blocksize)) && ((p->numentries) > 15))
450 p->issparse=0;
451 else
452 p->issparse=1;
453
454 // Note: everything is 1-indexed, so both locations and indices are shifted
455 for (int icol=1; icol < tmpBlock->startSize; icol++)
456 for (int jent=tmpBlock->start[icol-1]; jent<tmpBlock->start[icol]; jent++)
457 {
458 p->iindices[jent+1] = icol;
459 p->jindices[jent+1] = tmpBlock->index[jent] + 1;
460 p->entries [jent+1] = ((ConstantMatrixValues*)tmpBlock->value)->el[jent];
461 }
462
463 if (prev == NULL)
464 {
465 mconstraints[i].blocks = p;
466 }
467 else
468 {
469 prev->next = p;
470 prev->nextbyblock = p;
471 }
472 prev = p;
473 }
474 }
475 }
476
477 //garbage collection
478 if (blockOffset != NULL) delete [] blockOffset;
479 if (blockSize != NULL) delete [] blockSize;
480 if (mtxRef != NULL) delete [] mtxRef;
481 if (isdiag != NULL) delete [] isdiag;
482 if (mtxBlocks != NULL)
483 {
484 for (int i=0; i < ncon+1; i++)
485 {
486 if (mtxBlocks[i] != NULL) delete mtxBlocks[i];
487 mtxBlocks[i] = NULL;
488 }
489 delete []mtxBlocks;
490 }
491
492 this->bCallbuildSolverInstance = true;
493 return;
494 }
495
496 catch(const ErrorClass& eclass)
497 {
498 if (blockOffset != NULL) delete [] blockOffset;
499 if (blockSize != NULL) delete [] blockSize;
500 if (mtxRef != NULL) delete [] mtxRef;
501 if (isdiag != NULL) delete [] isdiag;
502 if (mtxBlocks != NULL)
503 {
504 for (int i=0; i <= ncon; i++)
505 {
506 if (mtxBlocks[i] != NULL) delete mtxBlocks[i];
507 mtxBlocks[i] = NULL;
508 }
509 delete []mtxBlocks;
510 }
511
512 osresult = new OSResult();
516 throw ErrorClass( osrl);
517 }
518
519}// end buildSolverInstance()
520
521
523{
534 struct paramstruc params;
535
536 std::ostringstream outStr;
537 std::ostringstream optStr;
538 int printlevel = 0;
539
540 try
541 {
542 /* get options from OSoL */
543 if(osoption == NULL && osol.length() > 0)
544 {
545 m_osolreader = new OSoLReader();
547 }
548
549 if( osoption != NULL && osoption->getNumberOfSolverOptions() > 0 )
550 {
551#ifndef NDEBUG
552 outStr.str("");
553 outStr.clear();
554 outStr << "number of solver options ";
556 outStr << std::endl;
558#endif
559
560 std::vector<SolverOption*> optionsVector;
561 optionsVector = osoption->getSolverOptions( "csdp",true);
562 char *pEnd;
563 int i;
564 int num_ipopt_options = optionsVector.size();
565 for(i = 0; i < num_ipopt_options; i++)
566 {
567#ifndef NDEBUG
568 outStr.str("");
569 outStr.clear();
570 outStr << "csdp solver option ";
571 outStr << optionsVector[ i]->name;
572 outStr << std::endl;
574#endif
575 optStr << optionsVector[ i]->name << "=" << optionsVector[ i]->value << std::endl;
576 }
577
578 FILE *paramfile;
579 paramfile=fopen("param.csdp","w");
580 if (!paramfile)
581 throw ErrorClass("File open error during option initialization");
582
583 fprintf(paramfile,"%s",(optStr.str()).c_str());
584 fclose(paramfile);
585 }
586 bSetSolverOptions = true;
587 return;
588 }
589 catch(const ErrorClass& eclass)
590 {
592 "Error in setSolverOption\n");
596 throw ErrorClass( osrl);
597 }
598}//setSolverOptions
599
600
601#if 0
602void CsdpSolver::setInitialValues()
603{
604 std::ostringstream outStr;
605 try
606 {
607 if(osinstance->getObjectiveNumber() <= 0)
608 throw ErrorClass("Ipopt NEEDS AN OBJECTIVE FUNCTION\n(For pure feasibility problems, use zero function.)");
609 this->bSetSolverOptions = true;
610 /* set the default options */
611 //app->Options()->SetNumericValue("tol", 1e-9);
612 app->Options()->SetIntegerValue("print_level", 0);
613 app->Options()->SetIntegerValue("max_iter", 20000);
614 app->Options()->SetNumericValue("bound_relax_factor", 0, true, true);
615 app->Options()->SetStringValue("mu_strategy", "adaptive", true, true);
616 //app->Options()->SetStringValue("output_file", "ipopt.out");
617 app->Options()->SetStringValue("check_derivatives_for_naninf", "yes");
618 // hessian constant for an LP
621 {
622 app->Options()->SetStringValue("hessian_constant", "yes", true, true);
623 }
625 {
626 if( osinstance->instanceData->objectives->obj[ 0]->maxOrMin.compare("min") != 0)
627 {
628 app->Options()->SetStringValue("nlp_scaling_method", "user-scaling");
629 }
630 }
631 /* end of the default options, now get options from OSoL */
632
633
634 if(osoption == NULL && osol.length() > 0)
635 {
636 m_osolreader = new OSoLReader();
638 }
639
640 if( osoption != NULL && osoption->getNumberOfSolverOptions() > 0 )
641 {
642#ifndef NDEBUG
643 outStr.str("");
644 outStr.clear();
645 outStr << "number of solver options ";
647 outStr << std::endl;
649#endif
650 std::vector<SolverOption*> optionsVector;
651 optionsVector = osoption->getSolverOptions( "ipopt",true);
652 char *pEnd;
653 int i;
654 int num_ipopt_options = optionsVector.size();
655 for(i = 0; i < num_ipopt_options; i++)
656 {
657#ifndef NDEBUG
658 outStr.str("");
659 outStr.clear();
660 outStr << "ipopt solver option ";
661 outStr << optionsVector[ i]->name;
662 outStr << std::endl;
664#endif
665 if(optionsVector[ i]->type == "numeric" )
666 {
667#ifndef NDEBUG
668 outStr.str("");
669 outStr.clear();
670 outStr << "FOUND A NUMERIC OPTION ";
671 outStr << os_strtod( optionsVector[ i]->value.c_str(), &pEnd );
672 outStr << std::endl;
674#endif
675 app->Options()->SetNumericValue(optionsVector[ i]->name, os_strtod( optionsVector[ i]->value.c_str(), &pEnd ) );
676 }
677 else if(optionsVector[ i]->type == "integer" )
678 {
679#ifndef NDEBUG
680 outStr.str("");
681 outStr.clear();
682 outStr << "FOUND AN INTEGER OPTION ";
683 outStr << atoi( optionsVector[ i]->value.c_str() );
684 outStr << std::endl;
686#endif
687 app->Options()->SetIntegerValue(optionsVector[ i]->name, atoi( optionsVector[ i]->value.c_str() ) );
688 }
689 else if(optionsVector[ i]->type == "string" )
690 {
691#ifndef NDEBUG
692 outStr.str("");
693 outStr.clear();
694 outStr << "FOUND A STRING OPTION ";
695 outStr << optionsVector[ i]->value.c_str();
696 outStr << std::endl;
698#endif
699 app->Options()->SetStringValue(optionsVector[ i]->name, optionsVector[ i]->value);
700 }
701 }
702 }
703 return;
704 }
705 catch(const ErrorClass& eclass)
706 {
707 osoutput->OSPrint(ENUM_OUTPUT_AREA_OSSolverInterfaces, ENUM_OUTPUT_LEVEL_error, "Error while setting options in IpoptSolver\n");
711 throw ErrorClass( osrl) ;
712 }
713}//end setInitialValues()
714
715#endif
716
717
719{
720 std::ostringstream outStr;
721
722 if( this->bCallbuildSolverInstance == false) buildSolverInstance();
723 if( this->bSetSolverOptions == false) setSolverOptions();
724 try
725 {
726
727// what about initial values for X and Y? perhaps even y?
728//if osoption->...->initialmatrix != NUll, set initial values. Make sure that defaults are there
729//in case X or Z is empty
730
731 /*
732 * Create an initial solution. This allocates space for X, y, and Z,
733 * and sets initial values.
734 */
735//else
737
738
739
740
741 //call solver
742 int returnCode = easy_sdp(nC_rows,ncon,C_matrix,rhsValues,mconstraints,0.0,&X,&y,&Z,&pobj,&dobj);
743
744 double* mdObjValues = NULL;
745 int solIdx = 0;
746 int numberOfOtherVariableResults;
747 int otherIdx;
748 int numCon = osinstance->getConstraintNumber();
749
751 {
752 mdObjValues = new double[1];
753 mdObjValues[0] = (dobj+pobj)/2;
754 outStr << std::endl << "Objective value f(x*) = " << os_dtoa_format(mdObjValues[0]);
756 }
757
758 std::string message = "Csdp solver finishes to the end.";
759 std::string solutionDescription = "";
760
761 // write resultHeader information
762 if(osresult->setSolverInvoked( "COIN-OR Csdp") != true)
763 throw ErrorClass("OSResult error: setSolverInvoked");
765 throw ErrorClass("OSResult error: setServiceName");
767 throw ErrorClass("OSResult error: setInstanceName");
768
769 //if(osresult->setJobID( osoption->jobID) != true)
770 // throw ErrorClass("OSResult error: setJobID");
771
772 // set basic problem parameters
774 throw ErrorClass("OSResult error: setVariableNumer");
775
776 if(osresult->setObjectiveNumber( 1) != true)
777 throw ErrorClass("OSResult error: setObjectiveNumber");
779 throw ErrorClass("OSResult error: setConstraintNumber");
780 if(osresult->setSolutionNumber( 1) != true)
781 throw ErrorClass("OSResult error: setSolutionNumer");
782 if(osresult->setGeneralMessage( message) != true)
783 throw ErrorClass("OSResult error: setGeneralMessage");
784
785
786 switch( returnCode)
787 {
788 case 0:
789 case 3:
790 {
791 if (returnCode == 0)
792 {
793 solutionDescription = "SUCCESS[Csdp]: Algorithm terminated normally at an optimal point, satisfying the convergence tolerances.";
794 osresult->setSolutionStatus(solIdx, "optimal", solutionDescription);
795 }
796 else
797 {
798 solutionDescription = "PARTIAL SUCCESS[Csdp]: A solution has been found, but full accuracy was not achieved.";
799 osresult->setSolutionStatus(solIdx, "unsure", solutionDescription);
800 }
801
803 osresult->setObjectiveValuesDense(solIdx, mdObjValues);
804 //set X, y, Z
805 if(numCon > 0)
806 osresult->setDualVariableValuesDense(solIdx, y+1); // !! Csdp uses Fortran indexing
807
808 int* colOffset = new int[X.nblocks+1];
809 colOffset[0] = 0;
810 for (int i=1; i<=X.nblocks; i++)
811 colOffset[i] = colOffset[i-1] + X.blocks[i].blocksize;
812
813 int* colOffsetD = new int[Z.nblocks+1];
814 colOffsetD[0] = 0;
815 for (int i=1; i<=Z.nblocks; i++)
816 colOffsetD[i] = colOffsetD[i-1] + Z.blocks[i].blocksize;
817
818 // initial the <variables> element: one primal matrixVar, one dual matrixVar (in <other>)
820 throw ErrorClass("OSResult error: setMatrixVariableSolution");
821
822
823 if (!osresult->setMatrixVarValuesAttributes(0,0,0,colOffset[X.nblocks],
824 colOffset[X.nblocks], ENUM_MATRIX_SYMMETRY_lower))
825 throw ErrorClass("OSResult error: setMatrixVarValuesAttributes");
826
828 "", "", "", "csdp", "", 1))
829 throw ErrorClass("OSResult error: setMatrixVariablesOtherResultGeneralAttributes");
831 colOffset[Z.nblocks],colOffset[Z.nblocks], ENUM_MATRIX_SYMMETRY_lower))
832 throw ErrorClass("OSResult error: setMatrixVariablesOtherResultMatrixAttributes");
833
834 if (!osresult->setMatrixVarValuesBlockStructure(0, 0, colOffset, X.nblocks + 1,
835 colOffset, X.nblocks + 1, X.nblocks))
836 throw ErrorClass("OSResult error: setMatrixVarValuesBlockStructure");
837
838 if (!osresult->setMatrixVariablesOtherResultBlockStructure(0, 0, 0, colOffsetD, Z.nblocks + 1,
839 colOffsetD, Z.nblocks + 1, Z.nblocks))
840 throw ErrorClass("OSResult error: setMatrixVariablesOtherResultBlockStructure");
841
842 int* start;
843 int* index;
845
846 int nonz;
847 double ent;
848
849 //count nonzeroes and set column starts
850 for (int blk=1; blk<=X.nblocks; blk++)
851 {
852 start = new int[colOffset[blk]-colOffset[blk-1]+1];
853 start[0] = 0;
854 nonz = 0;
855
856 switch (X.blocks[blk].blockcategory)
857 {
858 case DIAG:
859 for (int i=1; i<=X.blocks[blk].blocksize; i++)
860 {
861 ent=X.blocks[blk].data.vec[i];
862 if (ent != 0.0)
863 nonz++;
864 start[i] = nonz;
865 };
866 break;
867 case MATRIX:
868 for (int i=1; i<=X.blocks[blk].blocksize; i++)
869 {
870 for (int j=i; j<=X.blocks[blk].blocksize; j++)
871 {
872 ent=X.blocks[blk].data.mat[ijtok(i,j,X.blocks[blk].blocksize)];
873 if (ent != 0.0)
874 nonz++;
875 };
876 start[i] = nonz;
877 };
878 break;
879 case PACKEDMATRIX:
880 default:
881 throw ErrorClass("Invalid Block Type in CSDP solution");
882 }; // end switch
883
884 // go through nonzeros a second time and store values
885 index = new int[nonz];
886 value = new ConstantMatrixValues();
887 value->numberOfEl = nonz;
888 value->el = new double[nonz];
889
890 nonz = 0;
891 switch (X.blocks[blk].blockcategory)
892 {
893 case DIAG:
894 for (int i=1; i<=X.blocks[blk].blocksize; i++)
895 {
896 ent=X.blocks[blk].data.vec[i];
897 if (ent != 0.0)
898 {
899 index[nonz] = i-1;
900 value->el[nonz] = ent;
901 nonz++;
902 }
903 };
904 break;
905 case MATRIX:
906 for (int i=1; i<=X.blocks[blk].blocksize; i++)
907 {
908 for (int j=i; j<=X.blocks[blk].blocksize; j++)
909 {
910 ent=X.blocks[blk].data.mat[ijtok(i,j,X.blocks[blk].blocksize)];
911 if (ent != 0.0)
912 {
913 index[nonz] = j-1;
914 value->el[nonz] = ent;
915 nonz++;
916 }
917 };
918 start[i] = nonz;
919 };
920 break;
921 case PACKEDMATRIX:
922 default:
923 throw ErrorClass("Invalid Block Type in CSDP solution");
924 }; // end switch
925
926
927 if (!osresult->setMatrixVarValuesBlockElements(0, 0, blk-1, blk-1, blk-1, nonz,
928 start, index, value,
931 throw ErrorClass("OSResult error: setMatrixVarValuesBlockElements");
932 }; // end block
933
934 // now the dual variables
935 for (int blk=1; blk<=Z.nblocks; blk++)
936 {
937 start = new int[colOffsetD[blk]-colOffsetD[blk-1]+1];
938 start[0] = 0;
939 nonz = 0;
940
941 switch (Z.blocks[blk].blockcategory)
942 {
943 case DIAG:
944 for (int i=1; i<=Z.blocks[blk].blocksize; i++)
945 {
946 ent=Z.blocks[blk].data.vec[i];
947 if (ent != 0.0)
948 nonz++;
949 start[i] = nonz;
950 };
951
952 break;
953 case MATRIX:
954 for (int i=1; i<=Z.blocks[blk].blocksize; i++)
955 {
956 for (int j=i; j<=Z.blocks[blk].blocksize; j++)
957 {
958 ent=Z.blocks[blk].data.mat[ijtok(i,j,Z.blocks[blk].blocksize)];
959 if (ent != 0.0)
960 nonz++;
961 };
962 start[i] = nonz;
963 };
964 break;
965 case PACKEDMATRIX:
966 default:
967 throw ErrorClass("Invalid Block Type in CSDP solution");
968 }; // end switch
969
970 // go through nonzeros a second time and store values
971 index = new int[nonz];
972 value = new ConstantMatrixValues();
973 value->numberOfEl = nonz;
974 value->el = new double[nonz];
975
976 nonz = 0;
977 switch (Z.blocks[blk].blockcategory)
978 {
979 case DIAG:
980 for (int i=1; i<=Z.blocks[blk].blocksize; i++)
981 {
982 ent=Z.blocks[blk].data.vec[i];
983 if (ent != 0.0)
984 {
985 index[nonz] = i-1;
986 value->el[nonz] = ent;
987 nonz++;
988 }
989 };
990 break;
991 case MATRIX:
992 for (int i=1; i<=Z.blocks[blk].blocksize; i++)
993 {
994 for (int j=i; j<=Z.blocks[blk].blocksize; j++)
995 {
996 ent=Z.blocks[blk].data.mat[ijtok(i,j,Z.blocks[blk].blocksize)];
997 if (ent != 0.0)
998 {
999 index[nonz] = j-1;
1000 value->el[nonz] = ent;
1001 nonz++;
1002 }
1003 };
1004 start[i] = nonz;
1005 };
1006 break;
1007 case PACKEDMATRIX:
1008 default:
1009 throw ErrorClass("Invalid Block Type in CSDP solution");
1010 }; // end switch
1011
1012 if (!osresult->setMatrixVariablesOtherResultBlockElements(0, 0, 0, blk-1, blk-1, blk-1,
1013 nonz, start, index, value,
1016 throw ErrorClass("OSResult error: setMatrixVariablesOtherResultBlockElements");
1017 }; // end block
1018 break;
1019 }
1020 case 1:
1021 solutionDescription = "PRIMAL_INFEASIBILITY[Csdp]: Problem is primal infeasible.";
1022 osresult->setSolutionStatus(solIdx, "infeasible", solutionDescription);
1023 break;
1024
1025 case 2:
1026 solutionDescription = "DUAL_INFEASIBILITY[Csdp]: Problem is dual infeasible.";
1027 osresult->setSolutionStatus(solIdx, "infeasible", solutionDescription);
1028 break;
1029
1030 case 4:
1031 solutionDescription = "MAXITER_EXCEEDED[Csdp]: Maximum number of iterations exceeded.";
1032 osresult->setSolutionStatus(solIdx, "stoppedByLimit", solutionDescription);
1033 break;
1034
1035 case 5:
1036 solutionDescription = "STUCK AT EDGE[Csdp]: Stuck at edge of primal infeasibility.";
1037 osresult->setSolutionStatus(solIdx, "error", solutionDescription);
1038 break;
1039
1040 case 6:
1041 solutionDescription = "STUCK AT EDGE[Csdp]: Stuck at edge of dual infeasibility.";
1042 osresult->setSolutionStatus(solIdx, "error", solutionDescription);
1043 break;
1044
1045 case 7:
1046 solutionDescription = "LACK OF PROGRESS[Csdp]: Stopped due to lack of progress.";
1047 osresult->setSolutionStatus(solIdx, "error", solutionDescription);
1048 break;
1049
1050 case 8:
1051 solutionDescription = "SINGULARITY DETECTED[Csdp]: X, Z or O was singular.";
1052 osresult->setSolutionStatus(solIdx, "error", solutionDescription);
1053 break;
1054
1055 case 9:
1056 solutionDescription = "NaN or INF[Csdp]: Detected NaN or Infinity during computations.";
1057 osresult->setSolutionStatus(solIdx, "error", solutionDescription);
1058 break;
1059 }
1060
1061 if (returnCode != 0)
1062 throw ErrorClass("Csdp FAILED TO SOLVE THE PROBLEM");
1063
1065 }
1066
1067 catch (const ErrorClass& eclass)
1068 {
1069 outStr.str("");
1070 outStr.clear();
1071 outStr << "error in OSCsdpSolver routine solve():\n" << eclass.errormsg << endl;
1073
1075 osresult->setGeneralStatusType( "error");
1077 throw ErrorClass( osrl) ;
1078 }
1079}
1080
1082{
1083 int i0, j0, blk0;
1084 struct sparseblock *pp;
1085 std::ostringstream outStr;
1086
1087 outStr << std::endl << "Check problem data:" << std::endl << std::endl;
1088
1089 outStr << "Dimension of matrices (number of rows): n=" << nC_rows << std::endl;
1090 outStr << "Number of constraints (and matrices A_i): k=" << ncon << std::endl;
1091 for (i0=1; i0 <= ncon; i0++)
1092 {
1093 outStr << std::endl << "Right-hand side of constraint " << i0 << ": "
1094 << rhsValues[i0] << std::endl;
1095 outStr << std::endl << "Data for matrix A_" << i0 << ":" << std::endl;
1096 pp = mconstraints[i0].blocks;
1097 while (pp != NULL)
1098 {
1099 outStr << std::endl << "Block " << pp->blocknum << ":" << std::endl;;
1100 outStr << "Block size: " << pp->blocksize << std::endl;
1101 outStr << "Number of entries: " << pp->numentries << std::endl;
1102 for (j0=1; j0 <= pp->numentries; j0++)
1103 outStr << "Entry in row " << pp->iindices[j0] << ", col " << pp->jindices[j0]
1104 << " has value " << pp->entries[j0] << std::endl;
1105 pp = pp->next;
1106 }
1107 }
1108
1109 outStr << std::endl << "Data for matrix C:" << std::endl;
1110 outStr << "Number of blocks: " << C_matrix.nblocks << std::endl;
1111 for (blk0=1; blk0 <= C_matrix.nblocks; blk0++)
1112 {
1113 outStr << std::endl << "Data for block " << blk0 << ":" << std::endl;
1114 outStr << "Size: " << C_matrix.blocks[blk0].blocksize << std::endl;
1115 if (C_matrix.blocks[blk0].blockcategory == DIAG)
1116 {
1117 outStr << "Type: diagonal" << std::endl;
1118 for (i0=1; i0 <= C_matrix.blocks[blk0].blocksize; i0++)
1119 outStr << "Entry in row " << i0 << ", col " << i0 << " has value "
1120 << C_matrix.blocks[blk0].data.vec[i0] << std::endl;
1121 }
1122 else
1123 {
1124 outStr << "Type: dense" << std::endl;
1125 for (i0=1; i0 <= C_matrix.blocks[blk0].blocksize; i0++)
1126 for (j0=1; j0 <= C_matrix.blocks[blk0].blocksize; j0++)
1127 outStr << "Entry in row " << i0 << ", col " << j0 << " has value "
1128 << C_matrix.blocks[blk0].data.mat[ijtok(i0,j0,C_matrix.blocks[blk0].blocksize)]
1129 << std::endl;
1130 }
1131 }
1133}
1134
1135
const OSSmartPtr< OSOutput > osoutput
Definition OSOutput.cpp:39
std::string os_dtoa_format(double x)
std::string OSgetVersionInfo()
double os_strtod(const char *s00, char **se)
Definition OSdtoa.cpp:2541
to represent the nonzeros in a constantMatrix element
Definition OSMatrix.h:502
double pobj
struct blockmatrix X Z
double * y
OSiLReader * m_osilreader
m_osilreader is an OSiLReader object used to create an osinstance from an osil string if needed
double dobj
virtual void setSolverOptions()
The implementation of the virtual functions.
double * rhsValues
struct blockmatrix C_matrix
OSrLWriter * osrlwriter
virtual void solve()
solve results in an instance being read into the Csdp data structures and optimized
virtual ~CsdpSolver()
the CsdpSolver class destructor
virtual void buildSolverInstance()
The implementation of the virtual functions.
CsdpSolver()
the CsdpSolver class constructor
void dataEchoCheck()
use this for debugging, print out the instance that the solver thinks it has and compare this with th...
OSoLReader * m_osolreader
m_osolreader is an OSoLReader object used to create an osoption from an osol string if needed
struct constraintmatrix * mconstraints
std::string * csdpErrorMsg
std::string osol
osol holds the options for the solver
bool bSetSolverOptions
bSetSolverOptions is set to true if setSolverOptions has been called, false otherwise
std::string osrl
osrl holds the solution or result of the model
OSInstance * osinstance
osinstance holds the problem instance in-memory as an OSInstance object
bool bCallbuildSolverInstance
bCallbuildSolverInstance is set to true if buildSolverService has been called
std::string osil
osil holds the problem instance as a std::string
OSOption * osoption
osoption holds the solver options in-memory as an OSOption object
OSResult * osresult
osresult holds the solution or result of the model in-memory as an OSResult object
used for throwing exceptions.
std::string errormsg
errormsg is the error that is causing the exception to be thrown
a sparse matrix data structure for matrices that can hold nonconstant values and have block structure...
Definition OSMatrix.h:1769
GeneralSparseMatrix * getBlock(int rowIdx, int colIdx)
a method to retrieve a particular block from a collection
int inodeInt
inodeInt is the unique integer assigned to the OSnLNode or OSnLMNode in OSParameters....
Definition OSnLNode.h:62
OSnLMNode ** m_mMatrixChildren
m_mMatrixChildren holds all the matrix-valued operands, if any.
Definition OSnLNode.h:89
a sparse matrix data structure for matrices that can hold nonconstant values
Definition OSMatrix.h:1655
bool isDiagonal()
a method to determine whether the matrix is diagonal
MatrixElementValues * value
value holds a general array of value elements in the matrix, which could be constants,...
Definition OSMatrix.h:1723
int * index
index holds an integer array of rowIdx (or colIdx) elements in coefMatrix (AMatrix).
Definition OSMatrix.h:1708
int * start
start holds an integer array of start elements in the matrix, which points to the start of a column (...
Definition OSMatrix.h:1702
int valueSize
valueSize is the dimension of the index and value arrays
Definition OSMatrix.h:1696
int startSize
startSize is the dimension of the starts array
Definition OSMatrix.h:1691
Matrices * matrices
matrices is a pointer to a Matrices object
Objectives * objectives
objectives is a pointer to a Objectives object
OSMatrix ** matrix
matrix is a pointer to an array of OSMatrix object pointers
Definition OSInstance.h:499
int numberOfEl
each type of value is stored as an array named "el".
Definition OSMatrix.h:327
int numberOfRows
Definition OSMatrix.h:1904
int * getColumnPartition()
get the column partition of the matrix
int * getRowPartition()
get the row partition of the matrix
int getColumnPartitionSize()
get the size of the column partition of a matrix
ExpandedMatrixBlocks * getBlocks(int *rowPartition, int rowPartitionSize, int *colPartition, int colPartitionSize, bool rowMajor, bool appendToBlockArray)
A method to extract a block from a larger matrix The result is a sparse matrix object,...
int numberOfColumns
Definition OSMatrix.h:1905
int getRowPartitionSize()
get the size of the row partition of a matrix
double * getConstraintLowerBounds()
Get constraint lower bounds.
int getNumberOfQuadraticTerms()
Get the number of specified (usually nonzero) qTerms in the quadratic coefficients.
int getNumberOfMatrixVariables()
Get the number of matrix variables.
int getConstraintNumber()
Get number of constraints.
int getLinearConstraintCoefficientNumber()
Get number of specified (usually nonzero) linear constraint coefficient values.
int getMatrixNumber()
Get the number of matrices.
InstanceData * instanceData
A pointer to an InstanceData object.
int getNumberOfNonlinearExpressions()
Get number of nonlinear expressions.
int getVariableNumber()
Get number of variables.
std::string getInstanceName()
Get instance name.
ScalarExpressionTree * getNonlinearExpressionTree(int rowIdx)
Get the expression tree for a given row index.
char * getConstraintTypes()
Get constraint types.
std::string * getObjectiveMaxOrMins()
Get objective maxOrMins.
int getObjectiveNumber()
Get number of objectives.
a data structure to represent a matrix object (derived from MatrixType)
Definition OSMatrix.h:2186
virtual ENUM_MATRIX_TYPE getMatrixType()
std::vector< SolverOption * > getSolverOptions(std::string solver_name)
Get the options associated with a given solver.
int getNumberOfSolverOptions()
Get the number of solver options.
The Result Class.
Definition OSResult.h:2549
bool setGeneralMessage(std::string message)
Set the general message.
bool setSolutionNumber(int number)
set the number of solutions.
bool setMatrixVariablesOtherResultBlockElements(int solIdx, int otherIdx, int matrixVarIdx, int blkno, int blkRowIdx, int blkColIdx, int nz, int *start, int *index, MatrixElementValues *value, ENUM_MATRIX_TYPE valueType, ENUM_MATRIX_SYMMETRY symmetry=ENUM_MATRIX_SYMMETRY_none, bool rowMajor=false)
A method to set the elements within a block of a matrixVar associated with the [j]th "other" result i...
bool setInstanceName(std::string instanceName)
Set instance name.
bool setMatrixVariableSolution(int solIdx, int numberOfMatrixVar_, int numberOfOtherMatrixVariableResults_)
Set the [i]th optimization solution's MatrixVariableSolution, where i equals the given solution index...
bool setObjectiveValuesDense(int solIdx, double *objectiveValues)
Set the [i]th optimization solution's objective values, where i equals the given solution index.
bool setMatrixVarValuesAttributes(int solIdx, int idx, int matrixVarIdx, int numberOfRows, int numberOfColumns, ENUM_MATRIX_SYMMETRY symmetry=ENUM_MATRIX_SYMMETRY_none, ENUM_MATRIX_TYPE type=ENUM_MATRIX_TYPE_unknown, std::string name="")
A method to set general attributes for a matrixVar in the [i]th optimization solution,...
bool setMatrixVariablesOtherResultGeneralAttributes(int solIdx, int idx, std::string name, std::string description, std::string value, std::string type, std::string solver, std::string category, int numberOfMatrixVar=0, std::string matrixType="", int numberOfEnumerations=0, std::string enumType="")
A method to set general attributes for another (non-standard/solver specific) result associated with ...
bool setGeneralStatusType(std::string type)
Set the general status type, which can be: success, error, warning.
bool setObjectiveNumber(int objectiveNumber)
Set the objective number.
bool setMatrixVariablesOtherResultMatrixAttributes(int solIdx, int otherIdx, int matrixVarIdx, int numberOfRows, int numberOfColumns, ENUM_MATRIX_SYMMETRY symmetry=ENUM_MATRIX_SYMMETRY_none, ENUM_MATRIX_TYPE type=ENUM_MATRIX_TYPE_unknown, std::string name="")
A method to set attributes for a matrixVar in the [j]th other result associated with matrix variables...
bool setMatrixVariablesOtherResultBlockStructure(int solIdx, int otherIdx, int matrixVarIdx, int *colOffset, int colOffsetSize, int *rowOffset, int rowOffsetSize, int numberOfBlocks, int blocksConstructorIdx=0)
A method to set the block structure for the values of a matrixVar associated with the [j]th "other" r...
bool setServiceName(std::string serviceName)
Set service name.
bool setSolverInvoked(std::string solverInvoked)
Set solver invoked.
bool setVariableNumber(int variableNumber)
Set the variable number.
bool setDualVariableValuesDense(int solIdx, double *y)
Set the [i]th optimization solution's dual variable values, where i equals the given solution index.
bool setSolutionStatus(int solIdx, std::string type, std::string description)
Set the [i]th optimization solution status, where i equals the given solution index.
bool setConstraintNumber(int constraintNumber)
Set the constraint number.
bool setMatrixVarValuesBlockStructure(int solIdx, int idx, int *colOffset, int colOffsetSize, int *rowOffset, int rowOffsetSize, int numberOfBlocks, int blocksConstructorIdx=0)
A method to set the block structure for the values of a matrixVar in the [i]th optimization solution,...
bool setMatrixVarValuesBlockElements(int solIdx, int idx, int blkno, int blkRowIdx, int blkColIdx, int nz, int *start, int *index, MatrixElementValues *value, ENUM_MATRIX_TYPE valueType, ENUM_MATRIX_SYMMETRY symmetry=ENUM_MATRIX_SYMMETRY_none, bool rowMajor=false)
A method to set the elements within a block of a matrixVar in the [i]th optimization solution,...
Used to read an OSiL string.
Definition OSiLReader.h:38
OSInstance * readOSiL(const std::string &osil)
parse the OSiL model instance.
The OSnLMNode Class for nonlinear expressions involving matrices.
Definition OSnLNode.h:1761
The OSnLNode Class for nonlinear expressions.
Definition OSnLNode.h:180
Used to read an OSoL string.
Definition OSoLReader.h:38
OSOption * readOSoL(const std::string &osol)
parse the OSoL solver options.
Take an OSResult object and write a string that validates against OSrL.
Definition OSrLWriter.h:31
std::string writeOSrL(OSResult *theosresult)
create an osrl string from an OSResult object
std::string maxOrMin
declare the objective function to be a max or a min
Definition OSInstance.h:157
Objective ** obj
coef is pointer to an array of ObjCoef object pointers
Definition OSInstance.h:205
Used to hold part of the instance in memory.
OSnLNode * m_treeRoot
m_treeRoot holds the root node (of OSnLNode type) of the expression tree.
#define OS_MATRIX_TIMES
#define OS_MATRIX_REFERENCE
@ ENUM_OUTPUT_LEVEL_debug
@ ENUM_OUTPUT_LEVEL_trace
@ ENUM_OUTPUT_LEVEL_error
@ ENUM_OUTPUT_LEVEL_summary
@ ENUM_MATRIX_TYPE_constant
#define OS_MATRIX_VAR
@ ENUM_MATRIX_SYMMETRY_lower
#define OS_MATRIX_TRACE
@ ENUM_OUTPUT_AREA_OSSolverInterfaces