Context Navigation

source: trunk/Clp/src/ClpCholeskyWssmpKKT.cpp @ 1879

Last change on this file since 1879 was 1723, checked in by forrest, 9 years ago
out some printf statements
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 21.3 KB

Line
1	/* $Id: ClpCholeskyWssmpKKT.cpp 1723 2011-04-17 15:07:10Z forrest $ */
2	// Copyright (C) 2004, International Business Machines
3	// Corporation and others. All Rights Reserved.
4	// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7	#include "CoinPragma.hpp"
8	#include "CoinHelperFunctions.hpp"
9	#include "ClpHelperFunctions.hpp"
10
11	#include "ClpInterior.hpp"
12	#include "ClpCholeskyWssmpKKT.hpp"
13	#include "ClpQuadraticObjective.hpp"
14	#include "ClpMessage.hpp"
15
16	//#############################################################################
17	// Constructors / Destructor / Assignment
18	//#############################################################################
19
20	//-------------------------------------------------------------------
21	// Default Constructor
22	//-------------------------------------------------------------------
23	ClpCholeskyWssmpKKT::ClpCholeskyWssmpKKT (int denseThreshold)
24	: ClpCholeskyBase(denseThreshold)
25	{
26	type_ = 21;
27	}
28
29	//-------------------------------------------------------------------
30	// Copy constructor
31	//-------------------------------------------------------------------
32	ClpCholeskyWssmpKKT::ClpCholeskyWssmpKKT (const ClpCholeskyWssmpKKT & rhs)
33	: ClpCholeskyBase(rhs)
34	{
35	}
36
37
38	//-------------------------------------------------------------------
39	// Destructor
40	//-------------------------------------------------------------------
41	ClpCholeskyWssmpKKT::~ClpCholeskyWssmpKKT ()
42	{
43	}
44
45	//----------------------------------------------------------------
46	// Assignment operator
47	//-------------------------------------------------------------------
48	ClpCholeskyWssmpKKT &
49	ClpCholeskyWssmpKKT::operator=(const ClpCholeskyWssmpKKT& rhs)
50	{
51	if (this != &rhs) {
52	ClpCholeskyBase::operator=(rhs);
53	}
54	return *this;
55	}
56	//-------------------------------------------------------------------
57	// Clone
58	//-------------------------------------------------------------------
59	ClpCholeskyBase * ClpCholeskyWssmpKKT::clone() const
60	{
61	return new ClpCholeskyWssmpKKT(*this);
62	}
63	// At present I can't get wssmp to work as my libraries seem to be out of sync
64	// so I have linked in ekkwssmp which is an older version
65	#ifndef USE_EKKWSSMP
66	extern "C" {
67	void F77_FUNC(wsetmaxthrds,WSETMAXTHRDS)(const int* NTHREADS);
68
69	void F77_FUNC(wssmp,WSSMP)(const int* N, const int* IA,
70	const int* JA, const double* AVALS,
71	double* DIAG, int* PERM,
72	int* INVP, double* B,
73	const int* LDB, const int* NRHS,
74	double* AUX, const int* NAUX,
75	int* MRP, int* IPARM,
76	double* DPARM);
77	void F77_FUNC_(wsmp_clear,WSMP_CLEAR)(void);
78	}
79	#else
80	/* minimum needed for user */
81	typedef struct EKKModel EKKModel;
82	typedef struct EKKContext EKKContext;
83
84
85	extern "C" {
86	EKKContext * ekk_initializeContext();
87	void ekk_endContext(EKKContext * context);
88	EKKModel * ekk_newModel(EKKContext * env, const char * name);
89	int ekk_deleteModel(EKKModel * model);
90	}
91	static EKKModel * model = NULL;
92	static EKKContext * context = NULL;
93	extern "C" void ekkwssmp(EKKModel , int n,
94	int * columnStart , int * rowIndex , double * element,
95	double * diagonal , int * perm , int * invp ,
96	double * rhs , int * ldb , int * nrhs ,
97	double * aux , int * naux ,
98	int * mrp , int * iparm , double * dparm);
99	static void F77_FUNC(wssmp,WSSMP)( int n, int ia, int *ja,
100	double avals, double diag, int perm, int invp,
101	double b, int ldb, int nrhs, double aux, int *
102	naux, int mrp, int iparm, double *dparm)
103	{
104	if (!context) {
105	/* initialize OSL environment */
106	context = ekk_initializeContext();
107	model = ekk_newModel(context, "");
108	}
109	ekkwssmp(model, n, ia, ja,
110	avals, diag, perm, invp,
111	b, ldb, nrhs, aux,
112	naux, mrp, iparm, dparm);
113	//ekk_deleteModel(model);
114	//ekk_endContext(context);
115	}
116	#endif
117
118	/* Orders rows and saves pointer to model */
119	int
120	ClpCholeskyWssmpKKT::order(ClpInterior * model)
121	{
122	int numberRowsModel = model->numberRows();
123	int numberColumns = model->numberColumns();
124	int numberTotal = numberColumns + numberRowsModel;
125	numberRows_ = 2 * numberRowsModel + numberColumns;
126	rowsDropped_ = new char [numberRows_];
127	memset(rowsDropped_, 0, numberRows_);
128	numberRowsDropped_ = 0;
129	model_ = model;
130	CoinPackedMatrix * quadratic = NULL;
131	ClpQuadraticObjective * quadraticObj =
132	(dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
133	if (quadraticObj)
134	quadratic = quadraticObj->quadraticObjective();
135	int numberElements = model_->clpMatrix()->getNumElements();
136	numberElements = numberElements + 2 * numberRowsModel + numberTotal;
137	if (quadratic)
138	numberElements += quadratic->getNumElements();
139	// Space for starts
140	choleskyStart_ = new CoinBigIndex[numberRows_+1];
141	const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
142	const int * columnLength = model_->clpMatrix()->getVectorLengths();
143	const int * row = model_->clpMatrix()->getIndices();
144	//const double * element = model_->clpMatrix()->getElements();
145	// Now we have size - create arrays and fill in
146	try {
147	choleskyRow_ = new int [numberElements];
148	} catch (...) {
149	// no memory
150	delete [] choleskyStart_;
151	choleskyStart_ = NULL;
152	return -1;
153	}
154	try {
155	sparseFactor_ = new double[numberElements];
156	} catch (...) {
157	// no memory
158	delete [] choleskyRow_;
159	choleskyRow_ = NULL;
160	delete [] choleskyStart_;
161	choleskyStart_ = NULL;
162	return -1;
163	}
164	int iRow, iColumn;
165
166	sizeFactor_ = 0;
167	// matrix
168	if (!quadratic) {
169	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
170	choleskyStart_[iColumn] = sizeFactor_;
171	choleskyRow_[sizeFactor_++] = iColumn;
172	CoinBigIndex start = columnStart[iColumn];
173	CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
174	for (CoinBigIndex j = start; j < end; j++) {
175	choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
176	}
177	}
178	} else {
179	// Quadratic
180	const int * columnQuadratic = quadratic->getIndices();
181	const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
182	const int * columnQuadraticLength = quadratic->getVectorLengths();
183	//const double * quadraticElement = quadratic->getElements();
184	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
185	choleskyStart_[iColumn] = sizeFactor_;
186	choleskyRow_[sizeFactor_++] = iColumn;
187	for (CoinBigIndex j = columnQuadraticStart[iColumn];
188	j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
189	int jColumn = columnQuadratic[j];
190	if (jColumn > iColumn)
191	choleskyRow_[sizeFactor_++] = jColumn;
192	}
193	CoinBigIndex start = columnStart[iColumn];
194	CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
195	for (CoinBigIndex j = start; j < end; j++) {
196	choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
197	}
198	}
199	}
200	// slacks
201	for (; iColumn < numberTotal; iColumn++) {
202	choleskyStart_[iColumn] = sizeFactor_;
203	choleskyRow_[sizeFactor_++] = iColumn;
204	choleskyRow_[sizeFactor_++] = iColumn - numberColumns + numberTotal;
205	}
206	// Transpose - nonzero diagonal (may regularize)
207	for (iRow = 0; iRow < numberRowsModel; iRow++) {
208	choleskyStart_[iRow+numberTotal] = sizeFactor_;
209	// diagonal
210	choleskyRow_[sizeFactor_++] = iRow + numberTotal;
211	}
212	choleskyStart_[numberRows_] = sizeFactor_;
213	permuteInverse_ = new int [numberRows_];
214	permute_ = new int[numberRows_];
215	integerParameters_[0] = 0;
216	int i0 = 0;
217	int i1 = 1;
218	#ifndef USE_EKKWSSMP
219	int i2 = 1;
220	if (model->numberThreads() <= 0)
221	i2 = 1;
222	else
223	i2 = model->numberThreads();
224	F77_FUNC(wsetmaxthrds,WSETMAXTHRDS)(&i2);
225	#endif
226	F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
227	NULL, permute_, permuteInverse_, 0, &numberRows_, &i1,
228	NULL, &i0, NULL, integerParameters_, doubleParameters_);
229	integerParameters_[1] = 1; //order and symbolic
230	integerParameters_[2] = 2;
231	integerParameters_[3] = 0; //CSR
232	integerParameters_[4] = 0; //C style
233	integerParameters_[13] = 1; //reuse initial factorization space
234	integerParameters_[15+0] = 1; //ordering
235	integerParameters_[15+1] = 0;
236	integerParameters_[15+2] = 1;
237	integerParameters_[15+3] = 0;
238	integerParameters_[15+4] = 1;
239	doubleParameters_[10] = 1.0e-20;
240	doubleParameters_[11] = 1.0e-15;
241	#if 1
242	integerParameters_[1] = 2; //just symbolic
243	for (int iRow = 0; iRow < numberRows_; iRow++) {
244	permuteInverse_[iRow] = iRow;
245	permute_[iRow] = iRow;
246	}
247	#endif
248	F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
249	NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
250	NULL, &i0, NULL, integerParameters_, doubleParameters_);
251	//std::cout<<"Ordering and symbolic factorization took "<<doubleParameters_[0]<<std::endl;
252	if (integerParameters_[63]) {
253	std::cout << "wssmp returning error code of " << integerParameters_[63] << std::endl;
254	return 1;
255	}
256	std::cout << integerParameters_[23] << " elements in sparse Cholesky" << std::endl;
257	if (!integerParameters_[23]) {
258	for (int iRow = 0; iRow < numberRows_; iRow++) {
259	permuteInverse_[iRow] = iRow;
260	permute_[iRow] = iRow;
261	}
262	std::cout << "wssmp says no elements - fully dense? - switching to dense" << std::endl;
263	integerParameters_[1] = 2;
264	integerParameters_[2] = 2;
265	integerParameters_[7] = 1; // no permute
266	F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
267	NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
268	NULL, &i0, NULL, integerParameters_, doubleParameters_);
269	std::cout << integerParameters_[23] << " elements in dense Cholesky" << std::endl;
270	}
271	return 0;
272	}
273	/* Does Symbolic factorization given permutation.
274	This is called immediately after order. If user provides this then
275	user must provide factorize and solve. Otherwise the default factorization is used
276	returns non-zero if not enough memory */
277	int
278	ClpCholeskyWssmpKKT::symbolic()
279	{
280	return 0;
281	}
282	/* Factorize - filling in rowsDropped and returning number dropped */
283	int
284	ClpCholeskyWssmpKKT::factorize(const double * diagonal, int * rowsDropped)
285	{
286	int numberRowsModel = model_->numberRows();
287	int numberColumns = model_->numberColumns();
288	int numberTotal = numberColumns + numberRowsModel;
289	int newDropped = 0;
290	double largest = 0.0;
291	double smallest;
292	//perturbation
293	double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
294	perturbation = perturbation * perturbation;
295	if (perturbation > 1.0) {
296	#ifdef COIN_DEVELOP
297	//if (model_->model()->logLevel()&4)
298	std::cout << "large perturbation " << perturbation << std::endl;
299	#endif
300	perturbation = sqrt(perturbation);;
301	perturbation = 1.0;
302	}
303	// need to recreate every time
304	int iRow, iColumn;
305	const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
306	const int * columnLength = model_->clpMatrix()->getVectorLengths();
307	const int * row = model_->clpMatrix()->getIndices();
308	const double * element = model_->clpMatrix()->getElements();
309
310	CoinBigIndex numberElements = 0;
311	CoinPackedMatrix * quadratic = NULL;
312	ClpQuadraticObjective * quadraticObj =
313	(dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
314	if (quadraticObj)
315	quadratic = quadraticObj->quadraticObjective();
316	// matrix
317	if (!quadratic) {
318	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
319	choleskyStart_[iColumn] = numberElements;
320	double value = diagonal[iColumn];
321	if (fabs(value) > 1.0e-100) {
322	value = 1.0 / value;
323	largest = CoinMax(largest, fabs(value));
324	sparseFactor_[numberElements] = -value;
325	choleskyRow_[numberElements++] = iColumn;
326	CoinBigIndex start = columnStart[iColumn];
327	CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
328	for (CoinBigIndex j = start; j < end; j++) {
329	choleskyRow_[numberElements] = row[j] + numberTotal;
330	sparseFactor_[numberElements++] = element[j];
331	largest = CoinMax(largest, fabs(element[j]));
332	}
333	} else {
334	sparseFactor_[numberElements] = -1.0e100;
335	choleskyRow_[numberElements++] = iColumn;
336	}
337	}
338	} else {
339	// Quadratic
340	const int * columnQuadratic = quadratic->getIndices();
341	const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
342	const int * columnQuadraticLength = quadratic->getVectorLengths();
343	const double * quadraticElement = quadratic->getElements();
344	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
345	choleskyStart_[iColumn] = numberElements;
346	CoinBigIndex savePosition = numberElements;
347	choleskyRow_[numberElements++] = iColumn;
348	double value = diagonal[iColumn];
349	if (fabs(value) > 1.0e-100) {
350	value = 1.0 / value;
351	for (CoinBigIndex j = columnQuadraticStart[iColumn];
352	j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
353	int jColumn = columnQuadratic[j];
354	if (jColumn > iColumn) {
355	sparseFactor_[numberElements] = -quadraticElement[j];
356	choleskyRow_[numberElements++] = jColumn;
357	} else if (iColumn == jColumn) {
358	value += quadraticElement[j];
359	}
360	}
361	largest = CoinMax(largest, fabs(value));
362	sparseFactor_[savePosition] = -value;
363	CoinBigIndex start = columnStart[iColumn];
364	CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
365	for (CoinBigIndex j = start; j < end; j++) {
366	choleskyRow_[numberElements] = row[j] + numberTotal;
367	sparseFactor_[numberElements++] = element[j];
368	largest = CoinMax(largest, fabs(element[j]));
369	}
370	} else {
371	value = 1.0e100;
372	sparseFactor_[savePosition] = -value;
373	}
374	}
375	}
376	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
377	assert (sparseFactor_[choleskyStart_[iColumn]] < 0.0);
378	}
379	// slacks
380	for (iColumn = numberColumns; iColumn < numberTotal; iColumn++) {
381	choleskyStart_[iColumn] = numberElements;
382	double value = diagonal[iColumn];
383	if (fabs(value) > 1.0e-100) {
384	value = 1.0 / value;
385	largest = CoinMax(largest, fabs(value));
386	} else {
387	value = 1.0e100;
388	}
389	sparseFactor_[numberElements] = -value;
390	choleskyRow_[numberElements++] = iColumn;
391	choleskyRow_[numberElements] = iColumn - numberColumns + numberTotal;
392	sparseFactor_[numberElements++] = -1.0;
393	}
394	// Finish diagonal
395	double delta2 = model_->delta(); // add delta*delta to bottom
396	delta2 *= delta2;
397	for (iRow = 0; iRow < numberRowsModel; iRow++) {
398	choleskyStart_[iRow+numberTotal] = numberElements;
399	choleskyRow_[numberElements] = iRow + numberTotal;
400	sparseFactor_[numberElements++] = delta2;
401	}
402	choleskyStart_[numberRows_] = numberElements;
403	int i1 = 1;
404	int i0 = 0;
405	integerParameters_[1] = 3;
406	integerParameters_[2] = 3;
407	integerParameters_[10] = 2;
408	//integerParameters_[11]=1;
409	integerParameters_[12] = 2;
410	// LDLT
411	integerParameters_[30] = 1;
412	doubleParameters_[20] = 1.0e100;
413	double largest2 = largest * 1.0e-20;
414	largest = CoinMin(largest2, 1.0e-11);
415	doubleParameters_[10] = CoinMax(1.0e-20, largest);
416	if (doubleParameters_[10] > 1.0e-3)
417	integerParameters_[9] = 1;
418	else
419	integerParameters_[9] = 0;
420	#ifndef WSMP
421	// Set up LDL cutoff
422	integerParameters_[34] = numberTotal;
423	doubleParameters_[20] = 1.0e-15;
424	doubleParameters_[34] = 1.0e-12;
425	//printf("tol is %g\n",doubleParameters_[10]);
426	//doubleParameters_[10]=1.0e-17;
427	#endif
428	int * rowsDropped2 = new int[numberRows_];
429	CoinZeroN(rowsDropped2, numberRows_);
430	F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
431	NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
432	NULL, &i0, rowsDropped2, integerParameters_, doubleParameters_);
433	//std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
434	if (integerParameters_[9]) {
435	std::cout << "scaling applied" << std::endl;
436	}
437	newDropped = integerParameters_[20];
438	#if 1
439	// Should save adjustments in ..R_
440	int n1 = 0, n2 = 0;
441	double * primalR = model_->primalR();
442	double * dualR = model_->dualR();
443	for (iRow = 0; iRow < numberTotal; iRow++) {
444	if (rowsDropped2[iRow]) {
445	n1++;
446	//printf("row region1 %d dropped\n",iRow);
447	//rowsDropped_[iRow]=1;
448	rowsDropped_[iRow] = 0;
449	primalR[iRow] = doubleParameters_[20];
450	} else {
451	rowsDropped_[iRow] = 0;
452	primalR[iRow] = 0.0;
453	}
454	}
455	for (; iRow < numberRows_; iRow++) {
456	if (rowsDropped2[iRow]) {
457	n2++;
458	//printf("row region2 %d dropped\n",iRow);
459	//rowsDropped_[iRow]=1;
460	rowsDropped_[iRow] = 0;
461	dualR[iRow-numberTotal] = doubleParameters_[34];
462	} else {
463	rowsDropped_[iRow] = 0;
464	dualR[iRow-numberTotal] = 0.0;
465	}
466	}
467	//printf("%d rows dropped in region1, %d in region2\n",n1,n2);
468	#endif
469	delete [] rowsDropped2;
470	//if (integerParameters_[20])
471	//std::cout<<integerParameters_[20]<<" rows dropped"<<std::endl;
472	largest = doubleParameters_[3];
473	smallest = doubleParameters_[4];
474	if (model_->messageHandler()->logLevel() > 1)
475	std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
476	choleskyCondition_ = largest / smallest;
477	if (integerParameters_[63] < 0)
478	return -1; // out of memory
479	status_ = 0;
480	return 0;
481	}
482	/* Uses factorization to solve. */
483	void
484	ClpCholeskyWssmpKKT::solve (double * region)
485	{
486	abort();
487	}
488	/* Uses factorization to solve. */
489	void
490	ClpCholeskyWssmpKKT::solveKKT (double * region1, double * region2, const double * diagonal,
491	double diagonalScaleFactor)
492	{
493	int numberRowsModel = model_->numberRows();
494	int numberColumns = model_->numberColumns();
495	int numberTotal = numberColumns + numberRowsModel;
496	double * array = new double [numberRows_];
497	CoinMemcpyN(region1, numberTotal, array);
498	CoinMemcpyN(region2, numberRowsModel, array + numberTotal);
499	int i1 = 1;
500	int i0 = 0;
501	integerParameters_[1] = 4;
502	integerParameters_[2] = 4;
503	#if 0
504	integerParameters_[5] = 3;
505	doubleParameters_[5] = 1.0e-10;
506	integerParameters_[6] = 6;
507	#endif
508	F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
509	NULL, permute_, permuteInverse_, array, &numberRows_, &i1,
510	NULL, &i0, NULL, integerParameters_, doubleParameters_);
511	#if 0
512	int iRow;
513	for (iRow = 0; iRow < numberTotal; iRow++) {
514	if (rowsDropped_[iRow] && fabs(array[iRow]) > 1.0e-8) {
515	printf("row region1 %d dropped %g\n", iRow, array[iRow]);
516	}
517	}
518	for (; iRow < numberRows_; iRow++) {
519	if (rowsDropped_[iRow] && fabs(array[iRow]) > 1.0e-8) {
520	printf("row region2 %d dropped %g\n", iRow, array[iRow]);
521	}
522	}
523	#endif
524	CoinMemcpyN(array + numberTotal, numberRowsModel, region2);
525	#if 1
526	CoinMemcpyN(array, numberTotal, region1);
527	#else
528	multiplyAdd(region2, numberRowsModel, -1.0, array + numberColumns, 0.0);
529	CoinZeroN(array, numberColumns);
530	model_->clpMatrix()->transposeTimes(1.0, region2, array);
531	for (int iColumn = 0; iColumn < numberTotal; iColumn++)
532	region1[iColumn] = diagonal[iColumn] * (array[iColumn] - region1[iColumn]);
533	#endif
534	delete [] array;
535	#if 0
536	if (integerParameters_[5]) {
537	std::cout << integerParameters_[5] << " refinements ";
538	}
539	std::cout << doubleParameters_[6] << std::endl;
540	#endif
541	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: