Context Navigation

source: trunk/ClpCholeskyBase.cpp @ 629

Last change on this file since 629 was 629, checked in by forrest, 15 years ago
memory leaks
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 82.4 KB

Line
1	// Copyright (C) 2002, International Business Machines
2	// Corporation and others. All Rights Reserved.
3
4	#include "CoinPragma.hpp"
5
6	#include <iostream>
7
8	#include "ClpCholeskyBase.hpp"
9	#include "ClpInterior.hpp"
10	#include "ClpHelperFunctions.hpp"
11	#include "CoinHelperFunctions.hpp"
12	#include "CoinSort.hpp"
13	#include "ClpCholeskyDense.hpp"
14	#include "ClpMessage.hpp"
15	#include "ClpQuadraticObjective.hpp"
16
17	//#############################################################################
18	// Constructors / Destructor / Assignment
19	//#############################################################################
20
21	//-------------------------------------------------------------------
22	// Default Constructor
23	//-------------------------------------------------------------------
24	ClpCholeskyBase::ClpCholeskyBase (int denseThreshold) :
25	type_(0),
26	doKKT_(false),
27	goDense_(0.7),
28	choleskyCondition_(0.0),
29	model_(NULL),
30	numberTrials_(),
31	numberRows_(0),
32	status_(0),
33	rowsDropped_(NULL),
34	permuteInverse_(NULL),
35	permute_(NULL),
36	numberRowsDropped_(0),
37	sparseFactor_(NULL),
38	choleskyStart_(NULL),
39	choleskyRow_(NULL),
40	indexStart_(NULL),
41	diagonal_(NULL),
42	workDouble_(NULL),
43	link_(NULL),
44	workInteger_(NULL),
45	clique_(NULL),
46	sizeFactor_(0),
47	sizeIndex_(0),
48	firstDense_(0),
49	rowCopy_(NULL),
50	whichDense_(NULL),
51	denseColumn_(NULL),
52	dense_(NULL),
53	denseThreshold_(denseThreshold)
54	{
55	memset(integerParameters_,0,64*sizeof(int));
56	memset(doubleParameters_,0,64*sizeof(double));
57	}
58
59	//-------------------------------------------------------------------
60	// Copy constructor
61	//-------------------------------------------------------------------
62	ClpCholeskyBase::ClpCholeskyBase (const ClpCholeskyBase & rhs) :
63	type_(rhs.type_),
64	doKKT_(rhs.doKKT_),
65	goDense_(rhs.goDense_),
66	choleskyCondition_(rhs.choleskyCondition_),
67	model_(rhs.model_),
68	numberTrials_(rhs.numberTrials_),
69	numberRows_(rhs.numberRows_),
70	status_(rhs.status_),
71	numberRowsDropped_(rhs.numberRowsDropped_)
72	{
73	rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_,numberRows_);
74	permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_,numberRows_);
75	permute_ = ClpCopyOfArray(rhs.permute_,numberRows_);
76	sizeFactor_=rhs.sizeFactor_;
77	sizeIndex_ = rhs.sizeIndex_;
78	firstDense_ = rhs.firstDense_;
79	sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
80	choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
81	indexStart_ = ClpCopyOfArray(rhs.indexStart_,numberRows_);
82	choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,sizeIndex_);
83	diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_);
84	#if CLP_LONG_CHOLESKY!=1
85	workDouble_ = ClpCopyOfArray(rhs.workDouble_,numberRows_);
86	#else
87	// actually long double
88	workDouble_ = (double ) ClpCopyOfArray((longWork ) rhs.workDouble_,numberRows_);
89	#endif
90	link_ = ClpCopyOfArray(rhs.link_,numberRows_);
91	workInteger_ = ClpCopyOfArray(rhs.workInteger_,numberRows_);
92	clique_ = ClpCopyOfArray(rhs.clique_,numberRows_);
93	memcpy(integerParameters_,rhs.integerParameters_,64*sizeof(int));
94	memcpy(doubleParameters_,rhs.doubleParameters_,64*sizeof(double));
95	rowCopy_ = rhs.rowCopy_->clone();
96	whichDense_ = NULL;
97	denseColumn_=NULL;
98	dense_=NULL;
99	denseThreshold_ = rhs.denseThreshold_;
100	}
101
102	//-------------------------------------------------------------------
103	// Destructor
104	//-------------------------------------------------------------------
105	ClpCholeskyBase::~ClpCholeskyBase ()
106	{
107	delete [] rowsDropped_;
108	delete [] permuteInverse_;
109	delete [] permute_;
110	delete [] sparseFactor_;
111	delete [] choleskyStart_;
112	delete [] choleskyRow_;
113	delete [] indexStart_;
114	delete [] diagonal_;
115	delete [] workDouble_;
116	delete [] link_;
117	delete [] workInteger_;
118	delete [] clique_;
119	delete rowCopy_;
120	delete [] whichDense_;
121	delete [] denseColumn_;
122	delete dense_;
123	}
124
125	//----------------------------------------------------------------
126	// Assignment operator
127	//-------------------------------------------------------------------
128	ClpCholeskyBase &
129	ClpCholeskyBase::operator=(const ClpCholeskyBase& rhs)
130	{
131	if (this != &rhs) {
132	type_ = rhs.type_;
133	doKKT_ = rhs.doKKT_;
134	goDense_ = rhs.goDense_;
135	choleskyCondition_ = rhs.choleskyCondition_;
136	model_ = rhs.model_;
137	numberTrials_ = rhs.numberTrials_;
138	numberRows_ = rhs.numberRows_;
139	status_ = rhs.status_;
140	numberRowsDropped_ = rhs.numberRowsDropped_;
141	delete [] rowsDropped_;
142	delete [] permuteInverse_;
143	delete [] permute_;
144	delete [] sparseFactor_;
145	delete [] choleskyStart_;
146	delete [] choleskyRow_;
147	delete [] indexStart_;
148	delete [] diagonal_;
149	delete [] workDouble_;
150	delete [] link_;
151	delete [] workInteger_;
152	delete [] clique_;
153	delete rowCopy_;
154	delete [] whichDense_;
155	delete [] denseColumn_;
156	delete dense_;
157	rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_,numberRows_);
158	permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_,numberRows_);
159	permute_ = ClpCopyOfArray(rhs.permute_,numberRows_);
160	sizeFactor_=rhs.sizeFactor_;
161	sizeIndex_ = rhs.sizeIndex_;
162	firstDense_ = rhs.firstDense_;
163	sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
164	choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
165	choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
166	indexStart_ = ClpCopyOfArray(rhs.indexStart_,numberRows_);
167	choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,sizeIndex_);
168	diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_);
169	#if CLP_LONG_CHOLESKY!=1
170	workDouble_ = ClpCopyOfArray(rhs.workDouble_,numberRows_);
171	#else
172	// actually long double
173	workDouble_ = (double ) ClpCopyOfArray((longWork ) rhs.workDouble_,numberRows_);
174	#endif
175	link_ = ClpCopyOfArray(rhs.link_,numberRows_);
176	workInteger_ = ClpCopyOfArray(rhs.workInteger_,numberRows_);
177	clique_ = ClpCopyOfArray(rhs.clique_,numberRows_);
178	delete rowCopy_;
179	rowCopy_ = rhs.rowCopy_->clone();
180	whichDense_ = NULL;
181	denseColumn_=NULL;
182	dense_=NULL;
183	denseThreshold_ = rhs.denseThreshold_;
184	}
185	return *this;
186	}
187	// reset numberRowsDropped and rowsDropped.
188	void
189	ClpCholeskyBase::resetRowsDropped()
190	{
191	numberRowsDropped_=0;
192	memset(rowsDropped_,0,numberRows_);
193	}
194	/* Uses factorization to solve. - given as if KKT.
195	region1 is rows+columns, region2 is rows */
196	void
197	ClpCholeskyBase::solveKKT (double * region1, double * region2, const double * diagonal,
198	double diagonalScaleFactor)
199	{
200	if (!doKKT_) {
201	int iColumn;
202	int numberColumns = model_->numberColumns();
203	int numberTotal = numberRows_+numberColumns;
204	double * region1Save = new double[numberTotal];
205	for (iColumn=0;iColumn<numberTotal;iColumn++) {
206	region1[iColumn] *= diagonal[iColumn];
207	region1Save[iColumn]=region1[iColumn];
208	}
209	multiplyAdd(region1+numberColumns,numberRows_,-1.0,region2,1.0);
210	model_->clpMatrix()->times(1.0,region1,region2);
211	double maximumRHS = maximumAbsElement(region2,numberRows_);
212	double scale=1.0;
213	double unscale=1.0;
214	if (maximumRHS>1.0e-30) {
215	if (maximumRHS<=0.5) {
216	double factor=2.0;
217	while (maximumRHS<=0.5) {
218	maximumRHS*=factor;
219	scale*=factor;
220	} /* endwhile */
221	} else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
222	double factor=0.5;
223	while (maximumRHS>=2.0) {
224	maximumRHS*=factor;
225	scale*=factor;
226	} /* endwhile */
227	}
228	unscale=diagonalScaleFactor/scale;
229	} else {
230	//effectively zero
231	scale=0.0;
232	unscale=0.0;
233	}
234	multiplyAdd(NULL,numberRows_,0.0,region2,scale);
235	solve(region2);
236	multiplyAdd(NULL,numberRows_,0.0,region2,unscale);
237	multiplyAdd(region2,numberRows_,-1.0,region1+numberColumns,0.0);
238	CoinZeroN(region1,numberColumns);
239	model_->clpMatrix()->transposeTimes(1.0,region2,region1);
240	for (iColumn=0;iColumn<numberTotal;iColumn++)
241	region1[iColumn] = region1[iColumn]*diagonal[iColumn]-region1Save[iColumn];
242	delete [] region1Save;
243	} else {
244	// KKT
245	int numberRowsModel = model_->numberRows();
246	int numberColumns = model_->numberColumns();
247	int numberTotal = numberColumns + numberRowsModel;
248	double * array = new double [numberRows_];
249	CoinMemcpyN(region1,numberTotal,array);
250	CoinMemcpyN(region2,numberRowsModel,array+numberTotal);
251	solve(array);
252	int iRow;
253	for (iRow=0;iRow<numberTotal;iRow++) {
254	if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
255	printf("row region1 %d dropped %g\n",iRow,array[iRow]);
256	}
257	}
258	for (;iRow<numberRows_;iRow++) {
259	if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
260	printf("row region2 %d dropped %g\n",iRow,array[iRow]);
261	}
262	}
263	CoinMemcpyN(array+numberTotal,numberRowsModel,region2);
264	CoinMemcpyN(array,numberTotal,region1);
265	delete [] array;
266	}
267	}
268	//-------------------------------------------------------------------
269	// Clone
270	//-------------------------------------------------------------------
271	ClpCholeskyBase * ClpCholeskyBase::clone() const
272	{
273	return new ClpCholeskyBase(*this);
274	}
275	// Forms ADAT - returns nonzero if not enough memory
276	int
277	ClpCholeskyBase::preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT)
278	{
279	delete rowCopy_;
280	rowCopy_ = model_->clpMatrix()->reverseOrderedCopy();
281	if (!doKKT) {
282	numberRows_ = model_->numberRows();
283	rowsDropped_ = new char [numberRows_];
284	memset(rowsDropped_,0,numberRows_);
285	numberRowsDropped_=0;
286	// Space for starts
287	choleskyStart_ = new CoinBigIndex[numberRows_+1];
288	const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
289	const int * columnLength = model_->clpMatrix()->getVectorLengths();
290	const int * row = model_->clpMatrix()->getIndices();
291	const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
292	const int * rowLength = rowCopy_->getVectorLengths();
293	const int * column = rowCopy_->getIndices();
294	// We need two arrays for counts
295	int * which = new int [numberRows_];
296	int * used = new int[numberRows_+1];
297	CoinZeroN(used,numberRows_);
298	int iRow;
299	sizeFactor_=0;
300	int numberColumns = model_->numberColumns();
301	int numberDense=0;
302	//denseThreshold_=3;
303	if (denseThreshold_>0) {
304	delete [] whichDense_;
305	delete [] denseColumn_;
306	delete dense_;
307	whichDense_ = new char[numberColumns];
308	int iColumn;
309	used[numberRows_]=0;
310	for (iColumn=0;iColumn<numberColumns;iColumn++) {
311	int length = columnLength[iColumn];
312	used[length] += 1;
313	}
314	int nLong=0;
315	int stop = CoinMax(denseThreshold_/2,100);
316	for (iRow=numberRows_;iRow>=stop;iRow--) {
317	if (used[iRow])
318	printf("%d columns are of length %d\n",used[iRow],iRow);
319	nLong += used[iRow];
320	if (nLong>50\|\|nLong>(numberColumns>>2))
321	break;
322	}
323	CoinZeroN(used,numberRows_);
324	for (iColumn=0;iColumn<numberColumns;iColumn++) {
325	if (columnLength[iColumn]<denseThreshold_) {
326	whichDense_[iColumn]=0;
327	} else {
328	whichDense_[iColumn]=1;
329	numberDense++;
330	}
331	}
332	if (!numberDense\|\|numberDense>100) {
333	// free
334	delete [] whichDense_;
335	whichDense_=NULL;
336	denseColumn_=NULL;
337	dense_=NULL;
338	} else {
339	// space for dense columns
340	denseColumn_ = new longDouble [numberDense*numberRows_];
341	// dense cholesky
342	dense_ = new ClpCholeskyDense();
343	dense_->reserveSpace(NULL,numberDense);
344	printf("Taking %d columns as dense\n",numberDense);
345	}
346	}
347	int offset = includeDiagonal ? 0 : 1;
348	if (lowerTriangular)
349	offset=-offset;
350	for (iRow=0;iRow<numberRows_;iRow++) {
351	int number=0;
352	// make sure diagonal exists if includeDiagonal
353	if (!offset) {
354	which[0]=iRow;
355	used[iRow]=1;
356	number=1;
357	}
358	CoinBigIndex startRow=rowStart[iRow];
359	CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
360	if (lowerTriangular) {
361	for (CoinBigIndex k=startRow;k<endRow;k++) {
362	int iColumn=column[k];
363	if (!whichDense_\|\|!whichDense_[iColumn]) {
364	CoinBigIndex start=columnStart[iColumn];
365	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
366	for (CoinBigIndex j=start;j<end;j++) {
367	int jRow=row[j];
368	if (jRow<=iRow+offset) {
369	if (!used[jRow]) {
370	used[jRow]=1;
371	which[number++]=jRow;
372	}
373	}
374	}
375	}
376	}
377	} else {
378	for (CoinBigIndex k=startRow;k<endRow;k++) {
379	int iColumn=column[k];
380	if (!whichDense_\|\|!whichDense_[iColumn]) {
381	CoinBigIndex start=columnStart[iColumn];
382	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
383	for (CoinBigIndex j=start;j<end;j++) {
384	int jRow=row[j];
385	if (jRow>=iRow+offset) {
386	if (!used[jRow]) {
387	used[jRow]=1;
388	which[number++]=jRow;
389	}
390	}
391	}
392	}
393	}
394	}
395	sizeFactor_ += number;
396	int j;
397	for (j=0;j<number;j++)
398	used[which[j]]=0;
399	}
400	delete [] which;
401	// Now we have size - create arrays and fill in
402	try {
403	choleskyRow_ = new int [sizeFactor_];
404	}
405	catch (...) {
406	// no memory
407	delete [] choleskyStart_;
408	choleskyStart_=NULL;
409	return -1;
410	}
411	sizeFactor_=0;
412	which = choleskyRow_;
413	for (iRow=0;iRow<numberRows_;iRow++) {
414	int number=0;
415	// make sure diagonal exists if includeDiagonal
416	if (!offset) {
417	which[0]=iRow;
418	used[iRow]=1;
419	number=1;
420	}
421	choleskyStart_[iRow]=sizeFactor_;
422	CoinBigIndex startRow=rowStart[iRow];
423	CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
424	if (lowerTriangular) {
425	for (CoinBigIndex k=startRow;k<endRow;k++) {
426	int iColumn=column[k];
427	if (!whichDense_\|\|!whichDense_[iColumn]) {
428	CoinBigIndex start=columnStart[iColumn];
429	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
430	for (CoinBigIndex j=start;j<end;j++) {
431	int jRow=row[j];
432	if (jRow<=iRow+offset) {
433	if (!used[jRow]) {
434	used[jRow]=1;
435	which[number++]=jRow;
436	}
437	}
438	}
439	}
440	}
441	} else {
442	for (CoinBigIndex k=startRow;k<endRow;k++) {
443	int iColumn=column[k];
444	if (!whichDense_\|\|!whichDense_[iColumn]) {
445	CoinBigIndex start=columnStart[iColumn];
446	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
447	for (CoinBigIndex j=start;j<end;j++) {
448	int jRow=row[j];
449	if (jRow>=iRow+offset) {
450	if (!used[jRow]) {
451	used[jRow]=1;
452	which[number++]=jRow;
453	}
454	}
455	}
456	}
457	}
458	}
459	sizeFactor_ += number;
460	int j;
461	for (j=0;j<number;j++)
462	used[which[j]]=0;
463	// Sort
464	std::sort(which,which+number);
465	// move which on
466	which += number;
467	}
468	choleskyStart_[numberRows_]=sizeFactor_;
469	delete [] used;
470	return 0;
471	} else {
472	int numberRowsModel = model_->numberRows();
473	int numberColumns = model_->numberColumns();
474	int numberTotal = numberColumns + numberRowsModel;
475	numberRows_ = 2*numberRowsModel+numberColumns;
476	rowsDropped_ = new char [numberRows_];
477	memset(rowsDropped_,0,numberRows_);
478	numberRowsDropped_=0;
479	CoinPackedMatrix * quadratic = NULL;
480	ClpQuadraticObjective * quadraticObj =
481	(dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
482	if (quadraticObj)
483	quadratic = quadraticObj->quadraticObjective();
484	int numberElements = model_->clpMatrix()->getNumElements();
485	numberElements = numberElements+2*numberRowsModel+numberTotal;
486	if (quadratic)
487	numberElements += quadratic->getNumElements();
488	// Space for starts
489	choleskyStart_ = new CoinBigIndex[numberRows_+1];
490	const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
491	const int * columnLength = model_->clpMatrix()->getVectorLengths();
492	const int * row = model_->clpMatrix()->getIndices();
493	//const double * element = model_->clpMatrix()->getElements();
494	// Now we have size - create arrays and fill in
495	try {
496	choleskyRow_ = new int [numberElements];
497	}
498	catch (...) {
499	// no memory
500	delete [] choleskyStart_;
501	choleskyStart_=NULL;
502	return -1;
503	}
504	int iRow,iColumn;
505
506	sizeFactor_=0;
507	// matrix
508	if (lowerTriangular) {
509	if (!quadratic) {
510	for (iColumn=0;iColumn<numberColumns;iColumn++) {
511	choleskyStart_[iColumn]=sizeFactor_;
512	choleskyRow_[sizeFactor_++]=iColumn;
513	CoinBigIndex start=columnStart[iColumn];
514	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
515	if (!includeDiagonal)
516	start++;
517	for (CoinBigIndex j=start;j<end;j++) {
518	choleskyRow_[sizeFactor_++]=row[j]+numberTotal;
519	}
520	}
521	} else {
522	// Quadratic
523	const int * columnQuadratic = quadratic->getIndices();
524	const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
525	const int * columnQuadraticLength = quadratic->getVectorLengths();
526	//const double * quadraticElement = quadratic->getElements();
527	for (iColumn=0;iColumn<numberColumns;iColumn++) {
528	choleskyStart_[iColumn]=sizeFactor_;
529	if (includeDiagonal)
530	choleskyRow_[sizeFactor_++]=iColumn;
531	CoinBigIndex j;
532	for ( j=columnQuadraticStart[iColumn];
533	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
534	int jColumn = columnQuadratic[j];
535	if (jColumn>iColumn)
536	choleskyRow_[sizeFactor_++]=jColumn;
537	}
538	CoinBigIndex start=columnStart[iColumn];
539	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
540	for ( j=start;j<end;j++) {
541	choleskyRow_[sizeFactor_++]=row[j]+numberTotal;
542	}
543	}
544	}
545	// slacks
546	for (;iColumn<numberTotal;iColumn++) {
547	choleskyStart_[iColumn]=sizeFactor_;
548	if (includeDiagonal)
549	choleskyRow_[sizeFactor_++]=iColumn;
550	choleskyRow_[sizeFactor_++]=iColumn-numberColumns+numberTotal;
551	}
552	// Transpose - nonzero diagonal (may regularize)
553	for (iRow=0;iRow<numberRowsModel;iRow++) {
554	choleskyStart_[iRow+numberTotal]=sizeFactor_;
555	// diagonal
556	if (includeDiagonal)
557	choleskyRow_[sizeFactor_++]=iRow+numberTotal;
558	}
559	choleskyStart_[numberRows_]=sizeFactor_;
560	} else {
561	// transpose
562	ClpMatrixBase * rowCopy = model_->clpMatrix()->reverseOrderedCopy();
563	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
564	const int * rowLength = rowCopy->getVectorLengths();
565	const int * column = rowCopy->getIndices();
566	if (!quadratic) {
567	for (iColumn=0;iColumn<numberColumns;iColumn++) {
568	choleskyStart_[iColumn]=sizeFactor_;
569	if (includeDiagonal)
570	choleskyRow_[sizeFactor_++]=iColumn;
571	}
572	} else {
573	// Quadratic
574	// transpose
575	CoinPackedMatrix quadraticT;
576	quadraticT.reverseOrderedCopyOf(*quadratic);
577	const int * columnQuadratic = quadraticT.getIndices();
578	const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
579	const int * columnQuadraticLength = quadraticT.getVectorLengths();
580	//const double * quadraticElement = quadraticT.getElements();
581	for (iColumn=0;iColumn<numberColumns;iColumn++) {
582	choleskyStart_[iColumn]=sizeFactor_;
583	for (CoinBigIndex j=columnQuadraticStart[iColumn];
584	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
585	int jColumn = columnQuadratic[j];
586	if (jColumn<iColumn)
587	choleskyRow_[sizeFactor_++]=jColumn;
588	}
589	if (includeDiagonal)
590	choleskyRow_[sizeFactor_++]=iColumn;
591	}
592	}
593	int iRow;
594	// slacks
595	for (iRow=0;iRow<numberRowsModel;iRow++) {
596	choleskyStart_[iRow+numberColumns]=sizeFactor_;
597	if (includeDiagonal)
598	choleskyRow_[sizeFactor_++]=iRow+numberColumns;
599	}
600	for (iRow=0;iRow<numberRowsModel;iRow++) {
601	choleskyStart_[iRow+numberTotal]=sizeFactor_;
602	CoinBigIndex start=rowStart[iRow];
603	CoinBigIndex end=rowStart[iRow]+rowLength[iRow];
604	for (CoinBigIndex j=start;j<end;j++) {
605	choleskyRow_[sizeFactor_++]=column[j];
606	}
607	// slack
608	choleskyRow_[sizeFactor_++]=numberColumns+iRow;
609	if (includeDiagonal)
610	choleskyRow_[sizeFactor_++]=iRow+numberTotal;
611	}
612	choleskyStart_[numberRows_]=sizeFactor_;
613	}
614	}
615	return 0;
616	}
617	/* Orders rows and saves pointer to matrix.and model */
618	int
619	ClpCholeskyBase::order(ClpInterior * model)
620	{
621	model_=model;
622	int numberRowsModel = model_->numberRows();
623	int numberColumns = model_->numberColumns();
624	int numberTotal = numberColumns + numberRowsModel;
625	CoinPackedMatrix * quadratic = NULL;
626	ClpQuadraticObjective * quadraticObj =
627	(dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
628	if (quadraticObj)
629	quadratic = quadraticObj->quadraticObjective();
630	if (!doKKT_) {
631	numberRows_ = model->numberRows();
632	} else {
633	numberRows_ = 2*numberRowsModel+numberColumns;
634	}
635	rowsDropped_ = new char [numberRows_];
636	numberRowsDropped_=0;
637	rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
638	const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
639	const int * columnLength = model_->clpMatrix()->getVectorLengths();
640	const int * row = model_->clpMatrix()->getIndices();
641	const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
642	const int * rowLength = rowCopy_->getVectorLengths();
643	const int * column = rowCopy_->getIndices();
644	// We need arrays for counts
645	int * which = new int [numberRows_];
646	int * used = new int[numberRows_+1];
647	int * count = new int[numberRows_];
648	CoinZeroN(count,numberRows_);
649	CoinZeroN(used,numberRows_);
650	int iRow;
651	sizeFactor_=0;
652	permute_ = new int[numberRows_];
653	for (iRow=0;iRow<numberRows_;iRow++)
654	permute_[iRow]=iRow;
655	if (!doKKT_) {
656	int numberDense=0;
657	if (denseThreshold_>0) {
658	delete [] whichDense_;
659	delete [] denseColumn_;
660	delete dense_;
661	whichDense_ = new char[numberColumns];
662	int iColumn;
663	used[numberRows_]=0;
664	for (iColumn=0;iColumn<numberColumns;iColumn++) {
665	int length = columnLength[iColumn];
666	used[length] += 1;
667	}
668	int nLong=0;
669	int stop = CoinMax(denseThreshold_/2,100);
670	for (iRow=numberRows_;iRow>=stop;iRow--) {
671	if (used[iRow])
672	printf("%d columns are of length %d\n",used[iRow],iRow);
673	nLong += used[iRow];
674	if (nLong>50\|\|nLong>(numberColumns>>2))
675	break;
676	}
677	CoinZeroN(used,numberRows_);
678	for (iColumn=0;iColumn<numberColumns;iColumn++) {
679	if (columnLength[iColumn]<denseThreshold_) {
680	whichDense_[iColumn]=0;
681	} else {
682	whichDense_[iColumn]=1;
683	numberDense++;
684	}
685	}
686	if (!numberDense\|\|numberDense>100) {
687	// free
688	delete [] whichDense_;
689	whichDense_=NULL;
690	denseColumn_=NULL;
691	dense_=NULL;
692	} else {
693	// space for dense columns
694	denseColumn_ = new longDouble [numberDense*numberRows_];
695	// dense cholesky
696	dense_ = new ClpCholeskyDense();
697	dense_->reserveSpace(NULL,numberDense);
698	printf("Taking %d columns as dense\n",numberDense);
699	}
700	}
701	/*
702	Get row counts and size
703	*/
704	for (iRow=0;iRow<numberRows_;iRow++) {
705	int number=1;
706	// make sure diagonal exists
707	which[0]=iRow;
708	used[iRow]=1;
709	CoinBigIndex startRow=rowStart[iRow];
710	CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
711	for (CoinBigIndex k=startRow;k<endRow;k++) {
712	int iColumn=column[k];
713	if (!whichDense_\|\|!whichDense_[iColumn]) {
714	CoinBigIndex start=columnStart[iColumn];
715	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
716	for (CoinBigIndex j=start;j<end;j++) {
717	int jRow=row[j];
718	if (jRow<iRow) {
719	if (!used[jRow]) {
720	used[jRow]=1;
721	which[number++]=jRow;
722	count[jRow]++;
723	}
724	}
725	}
726	}
727	}
728	sizeFactor_ += number;
729	count[iRow]+=number;
730	int j;
731	for (j=0;j<number;j++)
732	used[which[j]]=0;
733	}
734	CoinSort_2(count,count+numberRows_,permute_);
735	} else {
736	// KKT
737	int numberElements = model_->clpMatrix()->getNumElements();
738	numberElements = numberElements+2*numberRowsModel+numberTotal;
739	if (quadratic)
740	numberElements += quadratic->getNumElements();
741	// off diagonal
742	numberElements -= numberRows_;
743	sizeFactor_=numberElements;
744	// If we sort we need to redo symbolic etc
745	}
746	delete [] which;
747	delete [] used;
748	delete [] count;
749	permuteInverse_ = new int [numberRows_];
750	memset(rowsDropped_,0,numberRows_);
751	for (iRow=0;iRow<numberRows_;iRow++) {
752	//permute_[iRow]=iRow; // force no permute
753	//permute_[iRow]=numberRows_-1-iRow; // force odd permute
754	//permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
755	permuteInverse_[permute_[iRow]]=iRow;
756	}
757	return 0;
758	}
759	/* Does Symbolic factorization given permutation.
760	This is called immediately after order. If user provides this then
761	user must provide factorize and solve. Otherwise the default factorization is used
762	returns non-zero if not enough memory */
763	int
764	ClpCholeskyBase::symbolic()
765	{
766	const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
767	const int * columnLength = model_->clpMatrix()->getVectorLengths();
768	const int * row = model_->clpMatrix()->getIndices();
769	const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
770	const int * rowLength = rowCopy_->getVectorLengths();
771	const int * column = rowCopy_->getIndices();
772	int numberRowsModel = model_->numberRows();
773	int numberColumns = model_->numberColumns();
774	int numberTotal = numberColumns + numberRowsModel;
775	CoinPackedMatrix * quadratic = NULL;
776	ClpQuadraticObjective * quadraticObj =
777	(dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
778	if (quadraticObj)
779	quadratic = quadraticObj->quadraticObjective();
780	// We need an array for counts
781	int * used = new int[numberRows_+1];
782	// If KKT then re-order so negative first
783	if (doKKT_) {
784	int nn=0;
785	int np=0;
786	int iRow;
787	for (iRow=0;iRow<numberRows_;iRow++) {
788	int originalRow = permute_[iRow];
789	if (originalRow<numberTotal)
790	permute_[nn++]=originalRow;
791	else
792	used[np++]=originalRow;
793	}
794	memcpy(permute_+nn,used,np*sizeof(int));
795	for (iRow=0;iRow<numberRows_;iRow++)
796	permuteInverse_[permute_[iRow]]=iRow;
797	}
798	CoinZeroN(used,numberRows_);
799	int iRow;
800	int iColumn;
801	bool noMemory=false;
802	CoinBigIndex * Astart = new CoinBigIndex[numberRows_+1];
803	int * Arow=NULL;
804	try {
805	Arow = new int [sizeFactor_];
806	}
807	catch (...) {
808	// no memory
809	delete [] Astart;
810	return -1;
811	}
812	choleskyStart_ = new int[numberRows_+1];
813	link_ = new int[numberRows_];
814	workInteger_ = new CoinBigIndex[numberRows_];
815	indexStart_ = new CoinBigIndex[numberRows_];
816	clique_ = new int[numberRows_];
817	// Redo so permuted upper triangular
818	sizeFactor_=0;
819	int * which = Arow;
820	if (!doKKT_) {
821	for (iRow=0;iRow<numberRows_;iRow++) {
822	int number=0;
823	int iOriginalRow = permute_[iRow];
824	Astart[iRow]=sizeFactor_;
825	CoinBigIndex startRow=rowStart[iOriginalRow];
826	CoinBigIndex endRow=rowStart[iOriginalRow]+rowLength[iOriginalRow];
827	for (CoinBigIndex k=startRow;k<endRow;k++) {
828	int iColumn=column[k];
829	if (!whichDense_\|\|!whichDense_[iColumn]) {
830	CoinBigIndex start=columnStart[iColumn];
831	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
832	for (CoinBigIndex j=start;j<end;j++) {
833	int jRow=row[j];
834	int jNewRow = permuteInverse_[jRow];
835	if (jNewRow<iRow) {
836	if (!used[jNewRow]) {
837	used[jNewRow]=1;
838	which[number++]=jNewRow;
839	}
840	}
841	}
842	}
843	}
844	sizeFactor_ += number;
845	int j;
846	for (j=0;j<number;j++)
847	used[which[j]]=0;
848	// Sort
849	std::sort(which,which+number);
850	// move which on
851	which += number;
852	}
853	} else {
854	// KKT
855	// transpose
856	ClpMatrixBase * rowCopy = model_->clpMatrix()->reverseOrderedCopy();
857	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
858	const int * rowLength = rowCopy->getVectorLengths();
859	const int * column = rowCopy->getIndices();
860	// temp
861	bool permuted=false;
862	for (iRow=0;iRow<numberRows_;iRow++) {
863	if (permute_[iRow]!=iRow) {
864	permuted=true;
865	break;
866	}
867	}
868	if (permuted) {
869	// Need to permute - ugly
870	if (!quadratic) {
871	for (iRow=0;iRow<numberRows_;iRow++) {
872	Astart[iRow]=sizeFactor_;
873	int iOriginalRow = permute_[iRow];
874	if (iOriginalRow<numberColumns) {
875	// A may be upper triangular by mistake
876	iColumn=iOriginalRow;
877	CoinBigIndex start=columnStart[iColumn];
878	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
879	for (CoinBigIndex j=start;j<end;j++) {
880	int kRow = row[j]+numberTotal;
881	kRow=permuteInverse_[kRow];
882	if (kRow<iRow)
883	Arow[sizeFactor_++]=kRow;
884	}
885	} else if (iOriginalRow<numberTotal) {
886	int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
887	if (kRow<iRow)
888	Arow[sizeFactor_++]=kRow;
889	} else {
890	int kRow = iOriginalRow-numberTotal;
891	CoinBigIndex start=rowStart[kRow];
892	CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
893	for (CoinBigIndex j=start;j<end;j++) {
894	int jRow = column[j];
895	int jNewRow = permuteInverse_[jRow];
896	if (jNewRow<iRow)
897	Arow[sizeFactor_++]=jNewRow;
898	}
899	// slack - should it be permute
900	kRow = permuteInverse_[kRow+numberColumns];
901	if (kRow<iRow)
902	Arow[sizeFactor_++]=kRow;
903	}
904	// Sort
905	std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
906	}
907	} else {
908	// quadratic
909	// transpose
910	CoinPackedMatrix quadraticT;
911	quadraticT.reverseOrderedCopyOf(*quadratic);
912	const int * columnQuadratic = quadraticT.getIndices();
913	const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
914	const int * columnQuadraticLength = quadraticT.getVectorLengths();
915	for (iRow=0;iRow<numberRows_;iRow++) {
916	Astart[iRow]=sizeFactor_;
917	int iOriginalRow = permute_[iRow];
918	if (iOriginalRow<numberColumns) {
919	// Quadratic bit
920	CoinBigIndex j;
921	for ( j=columnQuadraticStart[iOriginalRow];
922	j<columnQuadraticStart[iOriginalRow]+columnQuadraticLength[iOriginalRow];j++) {
923	int jColumn = columnQuadratic[j];
924	int jNewColumn = permuteInverse_[jColumn];
925	if (jNewColumn<iRow)
926	Arow[sizeFactor_++]=jNewColumn;
927	}
928	// A may be upper triangular by mistake
929	iColumn=iOriginalRow;
930	CoinBigIndex start=columnStart[iColumn];
931	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
932	for (j=start;j<end;j++) {
933	int kRow = row[j]+numberTotal;
934	kRow=permuteInverse_[kRow];
935	if (kRow<iRow)
936	Arow[sizeFactor_++]=kRow;
937	}
938	} else if (iOriginalRow<numberTotal) {
939	int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
940	if (kRow<iRow)
941	Arow[sizeFactor_++]=kRow;
942	} else {
943	int kRow = iOriginalRow-numberTotal;
944	CoinBigIndex start=rowStart[kRow];
945	CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
946	for (CoinBigIndex j=start;j<end;j++) {
947	int jRow = column[j];
948	int jNewRow = permuteInverse_[jRow];
949	if (jNewRow<iRow)
950	Arow[sizeFactor_++]=jNewRow;
951	}
952	// slack - should it be permute
953	kRow = permuteInverse_[kRow+numberColumns];
954	if (kRow<iRow)
955	Arow[sizeFactor_++]=kRow;
956	}
957	// Sort
958	std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
959	}
960	}
961	} else {
962	if (!quadratic) {
963	for (iRow=0;iRow<numberRows_;iRow++) {
964	Astart[iRow]=sizeFactor_;
965	}
966	} else {
967	// Quadratic
968	// transpose
969	CoinPackedMatrix quadraticT;
970	quadraticT.reverseOrderedCopyOf(*quadratic);
971	const int * columnQuadratic = quadraticT.getIndices();
972	const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
973	const int * columnQuadraticLength = quadraticT.getVectorLengths();
974	//const double * quadraticElement = quadraticT.getElements();
975	for (iRow=0;iRow<numberColumns;iRow++) {
976	int iOriginalRow = permute_[iRow];
977	Astart[iRow]=sizeFactor_;
978	for (CoinBigIndex j=columnQuadraticStart[iOriginalRow];
979	j<columnQuadraticStart[iOriginalRow]+columnQuadraticLength[iOriginalRow];j++) {
980	int jColumn = columnQuadratic[j];
981	int jNewColumn = permuteInverse_[jColumn];
982	if (jNewColumn<iRow)
983	Arow[sizeFactor_++]=jNewColumn;
984	}
985	}
986	}
987	int iRow;
988	// slacks
989	for (iRow=0;iRow<numberRowsModel;iRow++) {
990	Astart[iRow+numberColumns]=sizeFactor_;
991	}
992	for (iRow=0;iRow<numberRowsModel;iRow++) {
993	Astart[iRow+numberTotal]=sizeFactor_;
994	CoinBigIndex start=rowStart[iRow];
995	CoinBigIndex end=rowStart[iRow]+rowLength[iRow];
996	for (CoinBigIndex j=start;j<end;j++) {
997	Arow[sizeFactor_++]=column[j];
998	}
999	// slack
1000	Arow[sizeFactor_++]=numberColumns+iRow;
1001	}
1002	}
1003	delete rowCopy;
1004	}
1005	Astart[numberRows_]=sizeFactor_;
1006	firstDense_=numberRows_;
1007	symbolic1(Astart,Arow);
1008	// Now fill in indices
1009	try {
1010	// too big
1011	choleskyRow_ = new int[sizeFactor_];
1012	}
1013	catch (...) {
1014	// no memory
1015	noMemory=true;
1016	}
1017	double sizeFactor=sizeFactor_;
1018	if (!noMemory) {
1019	// Do lower triangular
1020	sizeFactor_=0;
1021	int * which = Arow;
1022	if (!doKKT_) {
1023	for (iRow=0;iRow<numberRows_;iRow++) {
1024	int number=0;
1025	int iOriginalRow = permute_[iRow];
1026	Astart[iRow]=sizeFactor_;
1027	if (!rowsDropped_[iOriginalRow]) {
1028	CoinBigIndex startRow=rowStart[iOriginalRow];
1029	CoinBigIndex endRow=rowStart[iOriginalRow]+rowLength[iOriginalRow];
1030	for (CoinBigIndex k=startRow;k<endRow;k++) {
1031	int iColumn=column[k];
1032	if (!whichDense_\|\|!whichDense_[iColumn]) {
1033	CoinBigIndex start=columnStart[iColumn];
1034	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1035	for (CoinBigIndex j=start;j<end;j++) {
1036	int jRow=row[j];
1037	int jNewRow = permuteInverse_[jRow];
1038	if (jNewRow>iRow&&!rowsDropped_[jRow]) {
1039	if (!used[jNewRow]) {
1040	used[jNewRow]=1;
1041	which[number++]=jNewRow;
1042	}
1043	}
1044	}
1045	}
1046	}
1047	sizeFactor_ += number;
1048	int j;
1049	for (j=0;j<number;j++)
1050	used[which[j]]=0;
1051	// Sort
1052	std::sort(which,which+number);
1053	// move which on
1054	which += number;
1055	}
1056	}
1057	} else {
1058	// KKT
1059	// temp
1060	bool permuted=false;
1061	for (iRow=0;iRow<numberRows_;iRow++) {
1062	if (permute_[iRow]!=iRow) {
1063	permuted=true;
1064	break;
1065	}
1066	}
1067	// but fake it
1068	for (iRow=0;iRow<numberRows_;iRow++) {
1069	//permute_[iRow]=iRow; // force no permute
1070	//permuteInverse_[permute_[iRow]]=iRow;
1071	}
1072	if (permuted) {
1073	// Need to permute - ugly
1074	if (!quadratic) {
1075	for (iRow=0;iRow<numberRows_;iRow++) {
1076	Astart[iRow]=sizeFactor_;
1077	int iOriginalRow = permute_[iRow];
1078	if (iOriginalRow<numberColumns) {
1079	iColumn=iOriginalRow;
1080	CoinBigIndex start=columnStart[iColumn];
1081	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1082	for (CoinBigIndex j=start;j<end;j++) {
1083	int kRow = row[j]+numberTotal;
1084	kRow=permuteInverse_[kRow];
1085	if (kRow>iRow)
1086	Arow[sizeFactor_++]=kRow;
1087	}
1088	} else if (iOriginalRow<numberTotal) {
1089	int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
1090	if (kRow>iRow)
1091	Arow[sizeFactor_++]=kRow;
1092	} else {
1093	int kRow = iOriginalRow-numberTotal;
1094	CoinBigIndex start=rowStart[kRow];
1095	CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
1096	for (CoinBigIndex j=start;j<end;j++) {
1097	int jRow = column[j];
1098	int jNewRow = permuteInverse_[jRow];
1099	if (jNewRow>iRow)
1100	Arow[sizeFactor_++]=jNewRow;
1101	}
1102	// slack - should it be permute
1103	kRow = permuteInverse_[kRow+numberColumns];
1104	if (kRow>iRow)
1105	Arow[sizeFactor_++]=kRow;
1106	}
1107	// Sort
1108	std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
1109	}
1110	} else {
1111	// quadratic
1112	const int * columnQuadratic = quadratic->getIndices();
1113	const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1114	const int * columnQuadraticLength = quadratic->getVectorLengths();
1115	for (iRow=0;iRow<numberRows_;iRow++) {
1116	Astart[iRow]=sizeFactor_;
1117	int iOriginalRow = permute_[iRow];
1118	if (iOriginalRow<numberColumns) {
1119	// Quadratic bit
1120	CoinBigIndex j;
1121	for ( j=columnQuadraticStart[iOriginalRow];
1122	j<columnQuadraticStart[iOriginalRow]+columnQuadraticLength[iOriginalRow];j++) {
1123	int jColumn = columnQuadratic[j];
1124	int jNewColumn = permuteInverse_[jColumn];
1125	if (jNewColumn>iRow)
1126	Arow[sizeFactor_++]=jNewColumn;
1127	}
1128	iColumn=iOriginalRow;
1129	CoinBigIndex start=columnStart[iColumn];
1130	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1131	for (j=start;j<end;j++) {
1132	int kRow = row[j]+numberTotal;
1133	kRow=permuteInverse_[kRow];
1134	if (kRow>iRow)
1135	Arow[sizeFactor_++]=kRow;
1136	}
1137	} else if (iOriginalRow<numberTotal) {
1138	int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
1139	if (kRow>iRow)
1140	Arow[sizeFactor_++]=kRow;
1141	} else {
1142	int kRow = iOriginalRow-numberTotal;
1143	CoinBigIndex start=rowStart[kRow];
1144	CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
1145	for (CoinBigIndex j=start;j<end;j++) {
1146	int jRow = column[j];
1147	int jNewRow = permuteInverse_[jRow];
1148	if (jNewRow>iRow)
1149	Arow[sizeFactor_++]=jNewRow;
1150	}
1151	// slack - should it be permute
1152	kRow = permuteInverse_[kRow+numberColumns];
1153	if (kRow>iRow)
1154	Arow[sizeFactor_++]=kRow;
1155	}
1156	// Sort
1157	std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
1158	}
1159	}
1160	} else {
1161	if (!quadratic) {
1162	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1163	Astart[iColumn]=sizeFactor_;
1164	CoinBigIndex start=columnStart[iColumn];
1165	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1166	for (CoinBigIndex j=start;j<end;j++) {
1167	Arow[sizeFactor_++]=row[j]+numberTotal;
1168	}
1169	}
1170	} else {
1171	// Quadratic
1172	const int * columnQuadratic = quadratic->getIndices();
1173	const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1174	const int * columnQuadraticLength = quadratic->getVectorLengths();
1175	//const double * quadraticElement = quadratic->getElements();
1176	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1177	Astart[iColumn]=sizeFactor_;
1178	CoinBigIndex j;
1179	for ( j=columnQuadraticStart[iColumn];
1180	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1181	int jColumn = columnQuadratic[j];
1182	if (jColumn>iColumn)
1183	Arow[sizeFactor_++]=jColumn;
1184	}
1185	CoinBigIndex start=columnStart[iColumn];
1186	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1187	for ( j=start;j<end;j++) {
1188	Arow[sizeFactor_++]=row[j]+numberTotal;
1189	}
1190	}
1191	}
1192	// slacks
1193	for (iRow=0;iRow<numberRowsModel;iRow++) {
1194	Astart[iRow+numberColumns]=sizeFactor_;
1195	Arow[sizeFactor_++]=iRow+numberTotal;
1196	}
1197	// Transpose - nonzero diagonal (may regularize)
1198	for (iRow=0;iRow<numberRowsModel;iRow++) {
1199	Astart[iRow+numberTotal]=sizeFactor_;
1200	}
1201	}
1202	Astart[numberRows_]=sizeFactor_;
1203	}
1204	symbolic2(Astart,Arow);
1205	if (sizeIndex_<sizeFactor_) {
1206	int * indices=NULL;
1207	try {
1208	indices = new int[sizeIndex_];
1209	}
1210	catch (...) {
1211	// no memory
1212	noMemory=true;
1213	}
1214	if (!noMemory) {
1215	memcpy(indices,choleskyRow_,sizeIndex_*sizeof(int));
1216	delete [] choleskyRow_;
1217	choleskyRow_=indices;
1218	}
1219	}
1220	}
1221	delete [] used;
1222	// Use cholesky regions
1223	delete [] Astart;
1224	delete [] Arow;
1225	double flops=0.0;
1226	for (iRow=0;iRow<numberRows_;iRow++) {
1227	int length = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1228	flops += ((double) length) * (length + 2.0);
1229	}
1230	if (model_->messageHandler()->logLevel()>0)
1231	std::cout<<sizeFactor<<" elements in sparse Cholesky, flop count "<<flops<<std::endl;
1232	try {
1233	sparseFactor_ = new longDouble [sizeFactor_];
1234	#if CLP_LONG_CHOLESKY!=1
1235	workDouble_ = new longDouble[numberRows_];
1236	#else
1237	// actually long double
1238	workDouble_ = (double *) (new longWork[numberRows_]);
1239	#endif
1240	diagonal_ = new longDouble[numberRows_];
1241	}
1242	catch (...) {
1243	// no memory
1244	noMemory=true;
1245	}
1246	if (noMemory) {
1247	delete [] choleskyRow_;
1248	choleskyRow_=NULL;
1249	delete [] choleskyStart_;
1250	choleskyStart_=NULL;
1251	delete [] permuteInverse_;
1252	permuteInverse_ = NULL;
1253	delete [] permute_;
1254	permute_ = NULL;
1255	delete [] choleskyStart_;
1256	choleskyStart_ = NULL;
1257	delete [] indexStart_;
1258	indexStart_ = NULL;
1259	delete [] link_;
1260	link_ = NULL;
1261	delete [] workInteger_;
1262	workInteger_ = NULL;
1263	delete [] sparseFactor_;
1264	sparseFactor_ = NULL;
1265	delete [] workDouble_;
1266	workDouble_ = NULL;
1267	delete [] diagonal_;
1268	diagonal_ = NULL;
1269	delete [] clique_;
1270	clique_ = NULL;
1271	return -1;
1272	}
1273	return 0;
1274	}
1275	int
1276	ClpCholeskyBase::symbolic1(const CoinBigIndex * Astart, const int * Arow)
1277	{
1278	int * marked = (int *) workInteger_;
1279	int iRow;
1280	// may not need to do this here but makes debugging easier
1281	for (iRow=0;iRow<numberRows_;iRow++) {
1282	marked[iRow]=-1;
1283	link_[iRow]=-1;
1284	choleskyStart_[iRow]=0; // counts
1285	}
1286	for (iRow=0;iRow<numberRows_;iRow++) {
1287	marked[iRow]=iRow;
1288	for (CoinBigIndex j=Astart[iRow];j<Astart[iRow+1];j++) {
1289	int kRow=Arow[j];
1290	while (marked[kRow] != iRow ) {
1291	if (link_[kRow] <0 )
1292	link_[kRow]=iRow;
1293	choleskyStart_[kRow]++;
1294	marked[kRow]=iRow;
1295	kRow = link_[kRow];
1296	}
1297	}
1298	}
1299	sizeFactor_=0;
1300	for (iRow=0;iRow<numberRows_;iRow++) {
1301	int number = choleskyStart_[iRow];
1302	choleskyStart_[iRow]=sizeFactor_;
1303	sizeFactor_ += number;
1304	}
1305	choleskyStart_[numberRows_]=sizeFactor_;
1306	return sizeFactor_;;
1307	}
1308	void
1309	ClpCholeskyBase::symbolic2(const CoinBigIndex * Astart, const int * Arow)
1310	{
1311	int * mergeLink = clique_;
1312	int * marker = (int *) workInteger_;
1313	int iRow;
1314	for (iRow=0;iRow<numberRows_;iRow++) {
1315	marker[iRow]=-1;
1316	mergeLink[iRow]=-1;
1317	link_[iRow]=-1; // not needed but makes debugging easier
1318	}
1319	int start=0;
1320	int end=0;
1321	choleskyStart_[0]=0;
1322
1323	for (iRow=0;iRow<numberRows_;iRow++) {
1324	int nz=0;
1325	int merge = mergeLink[iRow];
1326	bool marked=false;
1327	if (merge<0)
1328	marker[iRow]=iRow;
1329	else
1330	marker[iRow]=merge;
1331	start = end;
1332	int startSub=start;
1333	link_[iRow]=numberRows_;
1334	CoinBigIndex j;
1335	for ( j=Astart[iRow];j<Astart[iRow+1];j++) {
1336	int kRow=Arow[j];
1337	int k=iRow;
1338	int linked = link_[iRow];
1339	while (linked<=kRow) {
1340	k=linked;
1341	linked = link_[k];
1342	}
1343	nz++;
1344	link_[k]=kRow;
1345	link_[kRow]=linked;
1346	if (marker[kRow] != marker[iRow])
1347	marked=true;
1348	}
1349	bool reuse=false;
1350	// Check if we can re-use indices
1351	if (!marked && merge>=0&&mergeLink[merge]<0) {
1352	// can re-use all
1353	startSub = indexStart_[merge]+1;
1354	nz=choleskyStart_[merge+1]-(choleskyStart_[merge]+1);
1355	reuse=true;
1356	} else {
1357	// See if we can re-use any
1358	int k=mergeLink[iRow];
1359	int maxLength=0;
1360	while (k>=0) {
1361	int length = choleskyStart_[k+1]-(choleskyStart_[k]+1);
1362	int start = indexStart_[k]+1;
1363	int stop = start+length;
1364	if (length>maxLength) {
1365	maxLength = length;
1366	startSub = start;
1367	}
1368	int linked = iRow;
1369	for (CoinBigIndex j=start;j<stop;j++) {
1370	int kRow=choleskyRow_[j];
1371	int kk=linked;
1372	linked = link_[kk];
1373	while (linked<kRow) {
1374	kk=linked;
1375	linked = link_[kk];
1376	}
1377	if (linked!=kRow) {
1378	nz++;
1379	link_[kk]=kRow;
1380	link_[kRow]=linked;
1381	linked=kRow;
1382	}
1383	}
1384	k=mergeLink[k];
1385	}
1386	if (nz== maxLength)
1387	reuse=true; // can re-use
1388	}
1389	//reuse=false; //temp
1390	if (!reuse) {
1391	end += nz;
1392	startSub=start;
1393	int kRow = iRow;
1394	for (int j=start;j<end;j++) {
1395	kRow=link_[kRow];
1396	choleskyRow_[j]=kRow;
1397	assert (kRow<numberRows_);
1398	marker[kRow]=iRow;
1399	}
1400	marker[iRow]=iRow;
1401	}
1402	indexStart_[iRow]=startSub;
1403	choleskyStart_[iRow+1]=choleskyStart_[iRow]+nz;
1404	if (nz>1) {
1405	int kRow = choleskyRow_[startSub];
1406	mergeLink[iRow]=mergeLink[kRow];
1407	mergeLink[kRow]=iRow;
1408	}
1409	// should not be needed
1410	//std::sort(choleskyRow_+indexStart_[iRow]
1411	// ,choleskyRow_+indexStart_[iRow]+nz);
1412	#define CLP_DEBUG
1413	#ifdef CLP_DEBUG
1414	int last=-1;
1415	for ( j=indexStart_[iRow];j<indexStart_[iRow]+nz;j++) {
1416	int kRow=choleskyRow_[j];
1417	assert (kRow>last);
1418	last=kRow;
1419	}
1420	#endif
1421	}
1422	sizeFactor_ = choleskyStart_[numberRows_];
1423	sizeIndex_ = start;
1424	// find dense segment here
1425	int numberleft=numberRows_;
1426	for (iRow=0;iRow<numberRows_;iRow++) {
1427	CoinBigIndex left=sizeFactor_-choleskyStart_[iRow];
1428	double n=numberleft;
1429	double threshold = n(n-1.0)0.5*goDense_;
1430	if ((double) left >= threshold)
1431	break;
1432	numberleft--;
1433	}
1434	//iRow=numberRows_;
1435	int nDense = numberRows_-iRow;
1436	#define DENSE_THRESHOLD 8
1437	// don't do if dense columns
1438	if (nDense>=DENSE_THRESHOLD&&!dense_) {
1439	printf("Going dense for last %d rows\n",nDense);
1440	// make sure we don't disturb any indices
1441	CoinBigIndex k=0;
1442	for (int jRow=0;jRow<iRow;jRow++) {
1443	int nz=choleskyStart_[jRow+1]-choleskyStart_[jRow];
1444	k=CoinMax(k,indexStart_[jRow]+nz);
1445	}
1446	indexStart_[iRow]=k;
1447	int j;
1448	for (j=iRow+1;j<numberRows_;j++) {
1449	choleskyRow_[k++]=j;
1450	indexStart_[j]=k;
1451	}
1452	sizeIndex_=k;
1453	assert (k<=sizeFactor_); // can't happen with any reasonable defaults
1454	k=choleskyStart_[iRow];
1455	for (j=iRow+1;j<=numberRows_;j++) {
1456	k += numberRows_-j;
1457	choleskyStart_[j]=k;
1458	}
1459	// allow for blocked dense
1460	ClpCholeskyDense dense;
1461	sizeFactor_=choleskyStart_[iRow]+dense.space(nDense);
1462	firstDense_=iRow;
1463	if (doKKT_) {
1464	// redo permute so negative ones first
1465	int putN=firstDense_;
1466	int putP=0;
1467	int numberRowsModel = model_->numberRows();
1468	int numberColumns = model_->numberColumns();
1469	int numberTotal = numberColumns + numberRowsModel;
1470	for (iRow=firstDense_;iRow<numberRows_;iRow++) {
1471	int originalRow=permute_[iRow];
1472	if (originalRow<numberTotal)
1473	permute_[putN++]=originalRow;
1474	else
1475	permuteInverse_[putP++]=originalRow;
1476	}
1477	for (iRow=putN;iRow<numberRows_;iRow++) {
1478	permute_[iRow]=permuteInverse_[iRow-putN];
1479	}
1480	for (iRow=0;iRow<numberRows_;iRow++) {
1481	permuteInverse_[permute_[iRow]]=iRow;
1482	}
1483	}
1484	}
1485	// Clean up clique info
1486	for (iRow=0;iRow<numberRows_;iRow++)
1487	clique_[iRow]=0;
1488	int lastClique=-1;
1489	bool inClique=false;
1490	for (iRow=1;iRow<firstDense_;iRow++) {
1491	int sizeLast = choleskyStart_[iRow]-choleskyStart_[iRow-1];
1492	int sizeThis = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1493	if (indexStart_[iRow]==indexStart_[iRow-1]+1&&
1494	sizeThis==sizeLast-1&&
1495	sizeThis) {
1496	// in clique
1497	if (!inClique) {
1498	inClique=true;
1499	lastClique=iRow-1;
1500	}
1501	} else if (inClique) {
1502	int sizeClique=iRow-lastClique;
1503	for (int i=lastClique;i<iRow;i++) {
1504	clique_[i]=sizeClique;
1505	sizeClique--;
1506	}
1507	inClique=false;
1508	}
1509	}
1510	if (inClique) {
1511	int sizeClique=iRow-lastClique;
1512	for (int i=lastClique;i<iRow;i++) {
1513	clique_[i]=sizeClique;
1514	sizeClique--;
1515	}
1516	}
1517	//for (iRow=0;iRow<numberRows_;iRow++)
1518	//clique_[iRow]=0;
1519	}
1520	/* Factorize - filling in rowsDropped and returning number dropped */
1521	int
1522	ClpCholeskyBase::factorize(const double * diagonal, int * rowsDropped)
1523	{
1524	const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
1525	const int * columnLength = model_->clpMatrix()->getVectorLengths();
1526	const int * row = model_->clpMatrix()->getIndices();
1527	const double * element = model_->clpMatrix()->getElements();
1528	const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
1529	const int * rowLength = rowCopy_->getVectorLengths();
1530	const int * column = rowCopy_->getIndices();
1531	const double * elementByRow = rowCopy_->getElements();
1532	int numberColumns=model_->clpMatrix()->getNumCols();
1533	//perturbation
1534	longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
1535	//perturbation=perturbationperturbation100000000.0;
1536	if (perturbation>1.0) {
1537	//if (model_->model()->logLevel()&4)
1538	//std::cout <<"large perturbation "<<perturbation<<std::endl;
1539	perturbation=sqrt(perturbation);;
1540	perturbation=1.0;
1541	}
1542	int iRow;
1543	int iColumn;
1544	longDouble * work = workDouble_;
1545	CoinZeroN(work,numberRows_);
1546	int newDropped=0;
1547	double largest=1.0;
1548	double smallest=COIN_DBL_MAX;
1549	int numberDense=0;
1550	if (!doKKT_) {
1551	const double * diagonalSlack = diagonal + numberColumns;
1552	if (dense_)
1553	numberDense=dense_->numberRows();
1554	if (whichDense_) {
1555	longDouble * denseDiagonal = dense_->diagonal();
1556	longDouble * dense = denseColumn_;
1557	int iDense=0;
1558	for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1559	if (whichDense_[iColumn]) {
1560	CoinZeroN(dense,numberRows_);
1561	CoinBigIndex start=columnStart[iColumn];
1562	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1563	if (diagonal[iColumn]) {
1564	denseDiagonal[iDense++]=1.0/diagonal[iColumn];
1565	for (CoinBigIndex j=start;j<end;j++) {
1566	int jRow=row[j];
1567	dense[jRow] = element[j];
1568	}
1569	} else {
1570	denseDiagonal[iDense++]=1.0;
1571	}
1572	dense += numberRows_;
1573	}
1574	}
1575	}
1576	longDouble delta2 = model_->delta(); // add delta*delta to diagonal
1577	delta2 *= delta2;
1578	// largest in initial matrix
1579	double largest2=1.0e-20;
1580	for (iRow=0;iRow<numberRows_;iRow++) {
1581	longDouble * put = sparseFactor_+choleskyStart_[iRow];
1582	int * which = choleskyRow_+indexStart_[iRow];
1583	int iOriginalRow = permute_[iRow];
1584	int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1585	if (!rowLength[iOriginalRow])
1586	rowsDropped_[iOriginalRow]=1;
1587	if (!rowsDropped_[iOriginalRow]) {
1588	CoinBigIndex startRow=rowStart[iOriginalRow];
1589	CoinBigIndex endRow=rowStart[iOriginalRow]+rowLength[iOriginalRow];
1590	work[iRow] = diagonalSlack[iOriginalRow]+delta2;
1591	for (CoinBigIndex k=startRow;k<endRow;k++) {
1592	int iColumn=column[k];
1593	if (!whichDense_\|\|!whichDense_[iColumn]) {
1594	CoinBigIndex start=columnStart[iColumn];
1595	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1596	longDouble multiplier = diagonal[iColumn]*elementByRow[k];
1597	for (CoinBigIndex j=start;j<end;j++) {
1598	int jRow=row[j];
1599	int jNewRow = permuteInverse_[jRow];
1600	if (jNewRow>=iRow&&!rowsDropped_[jRow]) {
1601	longDouble value=element[j]*multiplier;
1602	work[jNewRow] += value;
1603	}
1604	}
1605	}
1606	}
1607	diagonal_[iRow]=work[iRow];
1608	largest2 = CoinMax(largest2,fabs(work[iRow]));
1609	work[iRow]=0.0;
1610	int j;
1611	for (j=0;j<number;j++) {
1612	int jRow = which[j];
1613	put[j]=work[jRow];
1614	largest2 = CoinMax(largest2,fabs(work[jRow]));
1615	work[jRow]=0.0;
1616	}
1617	} else {
1618	// dropped
1619	diagonal_[iRow]=1.0;
1620	int j;
1621	for (j=1;j<number;j++) {
1622	put[j]=0.0;
1623	}
1624	}
1625	}
1626	//check sizes
1627	largest2*=1.0e-20;
1628	largest = CoinMin(largest2,1.0e-11);
1629	int numberDroppedBefore=0;
1630	for (iRow=0;iRow<numberRows_;iRow++) {
1631	int dropped=rowsDropped_[iRow];
1632	// Move to int array
1633	rowsDropped[iRow]=dropped;
1634	if (!dropped) {
1635	longDouble diagonal = diagonal_[iRow];
1636	if (diagonal>largest2) {
1637	diagonal_[iRow]=diagonal+perturbation;
1638	} else {
1639	diagonal_[iRow]=diagonal+perturbation;
1640	rowsDropped[iRow]=2;
1641	numberDroppedBefore++;
1642	//printf("dropped - small diagonal %g\n",diagonal);
1643	}
1644	}
1645	}
1646	doubleParameters_[10]=CoinMax(1.0e-20,largest);
1647	integerParameters_[20]=0;
1648	doubleParameters_[3]=0.0;
1649	doubleParameters_[4]=COIN_DBL_MAX;
1650	integerParameters_[34]=0; // say all must be positive
1651	factorizePart2(rowsDropped);
1652	newDropped=integerParameters_[20]+numberDroppedBefore;
1653	largest=doubleParameters_[3];
1654	smallest=doubleParameters_[4];
1655	if (model_->messageHandler()->logLevel()>1)
1656	std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
1657	choleskyCondition_=largest/smallest;
1658	if (whichDense_) {
1659	int i;
1660	for ( i=0;i<numberRows_;i++) {
1661	assert (diagonal_[i]>=0.0);
1662	diagonal_[i]=sqrt(diagonal_[i]);
1663	}
1664	// Update dense columns (just L)
1665	// Zero out dropped rows
1666	for (i=0;i<numberDense;i++) {
1667	longDouble * a = denseColumn_+i*numberRows_;
1668	for (int j=0;j<numberRows_;j++) {
1669	if (rowsDropped[j])
1670	a[j]=0.0;
1671	}
1672	solveLong(a,1);
1673	}
1674	dense_->resetRowsDropped();
1675	longDouble * denseBlob = dense_->aMatrix();
1676	longDouble * denseDiagonal = dense_->diagonal();
1677	// Update dense matrix
1678	for (i=0;i<numberDense;i++) {
1679	const longDouble * a = denseColumn_+i*numberRows_;
1680	// do diagonal
1681	longDouble value = denseDiagonal[i];
1682	const longDouble * b = denseColumn_+i*numberRows_;
1683	for (int k=0;k<numberRows_;k++)
1684	value += a[k]*b[k];
1685	denseDiagonal[i]=value;
1686	for (int j=i+1;j<numberDense;j++) {
1687	longDouble value = 0.0;
1688	const longDouble * b = denseColumn_+j*numberRows_;
1689	for (int k=0;k<numberRows_;k++)
1690	value += a[k]*b[k];
1691	*denseBlob=value;
1692	denseBlob++;
1693	}
1694	}
1695	// dense cholesky (? long double)
1696	int * dropped = new int [numberDense];
1697	dense_->factorizePart2(dropped);
1698	delete [] dropped;
1699	}
1700	// try allowing all every time
1701	//printf("trying ?\n");
1702	//for (iRow=0;iRow<numberRows_;iRow++) {
1703	//rowsDropped[iRow]=0;
1704	//rowsDropped_[iRow]=0;
1705	//}
1706	bool cleanCholesky;
1707	//if (model_->numberIterations()<20\|\|(model_->numberIterations()&1)==0)
1708	if (model_->numberIterations()<2000)
1709	cleanCholesky=true;
1710	else
1711	cleanCholesky=false;
1712	if (cleanCholesky) {
1713	//drop fresh makes some formADAT easier
1714	if (newDropped\|\|numberRowsDropped_) {
1715	newDropped=0;
1716	for (int i=0;i<numberRows_;i++) {
1717	char dropped = rowsDropped[i];
1718	rowsDropped_[i]=dropped;
1719	rowsDropped_[i]=0;
1720	if (dropped==2) {
1721	//dropped this time
1722	rowsDropped[newDropped++]=i;
1723	rowsDropped_[i]=0;
1724	}
1725	}
1726	numberRowsDropped_=newDropped;
1727	newDropped=-(2+newDropped);
1728	}
1729	} else {
1730	if (newDropped) {
1731	newDropped=0;
1732	for (int i=0;i<numberRows_;i++) {
1733	char dropped = rowsDropped[i];
1734	rowsDropped_[i]=dropped;
1735	if (dropped==2) {
1736	//dropped this time
1737	rowsDropped[newDropped++]=i;
1738	rowsDropped_[i]=1;
1739	}
1740	}
1741	}
1742	numberRowsDropped_+=newDropped;
1743	if (numberRowsDropped_&&0) {
1744	std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
1745	numberRowsDropped_<<" dropped)";
1746	if (newDropped) {
1747	std::cout<<" ( "<<newDropped<<" dropped this time)";
1748	}
1749	std::cout<<std::endl;
1750	}
1751	}
1752	} else {
1753	//KKT
1754	CoinPackedMatrix * quadratic = NULL;
1755	ClpQuadraticObjective * quadraticObj =
1756	(dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
1757	if (quadraticObj)
1758	quadratic = quadraticObj->quadraticObjective();
1759	// matrix
1760	int numberRowsModel = model_->numberRows();
1761	int numberColumns = model_->numberColumns();
1762	int numberTotal = numberColumns + numberRowsModel;
1763	// temp
1764	bool permuted=false;
1765	for (iRow=0;iRow<numberRows_;iRow++) {
1766	if (permute_[iRow]!=iRow) {
1767	permuted=true;
1768	break;
1769	}
1770	}
1771	// but fake it
1772	for (iRow=0;iRow<numberRows_;iRow++) {
1773	//permute_[iRow]=iRow; // force no permute
1774	//permuteInverse_[permute_[iRow]]=iRow;
1775	}
1776	if (permuted) {
1777	longDouble delta2 = model_->delta(); // add delta*delta to bottom
1778	delta2 *= delta2;
1779	// Need to permute - ugly
1780	if (!quadratic) {
1781	for (iRow=0;iRow<numberRows_;iRow++) {
1782	longDouble * put = sparseFactor_+choleskyStart_[iRow];
1783	int * which = choleskyRow_+indexStart_[iRow];
1784	int iOriginalRow = permute_[iRow];
1785	if (iOriginalRow<numberColumns) {
1786	iColumn=iOriginalRow;
1787	longDouble value = diagonal[iColumn];
1788	if (fabs(value)>1.0e-100) {
1789	value = 1.0/value;
1790	largest = CoinMax(largest,fabs(value));
1791	diagonal_[iRow] = -value;
1792	CoinBigIndex start=columnStart[iColumn];
1793	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1794	for (CoinBigIndex j=start;j<end;j++) {
1795	int kRow = row[j]+numberTotal;
1796	kRow=permuteInverse_[kRow];
1797	if (kRow>iRow) {
1798	work[kRow]=element[j];
1799	largest = CoinMax(largest,fabs(element[j]));
1800	}
1801	}
1802	} else {
1803	diagonal_[iRow] = -value;
1804	}
1805	} else if (iOriginalRow<numberTotal) {
1806	longDouble value = diagonal[iOriginalRow];
1807	if (fabs(value)>1.0e-100) {
1808	value = 1.0/value;
1809	largest = CoinMax(largest,fabs(value));
1810	} else {
1811	value = 1.0e100;
1812	}
1813	diagonal_[iRow] = -value;
1814	int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
1815	if (kRow>iRow)
1816	work[kRow]=-1.0;
1817	} else {
1818	diagonal_[iRow]=delta2;
1819	int kRow = iOriginalRow-numberTotal;
1820	CoinBigIndex start=rowStart[kRow];
1821	CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
1822	for (CoinBigIndex j=start;j<end;j++) {
1823	int jRow = column[j];
1824	int jNewRow = permuteInverse_[jRow];
1825	if (jNewRow>iRow) {
1826	work[jNewRow]=elementByRow[j];
1827	largest = CoinMax(largest,fabs(elementByRow[j]));
1828	}
1829	}
1830	// slack - should it be permute
1831	kRow = permuteInverse_[kRow+numberColumns];
1832	if (kRow>iRow)
1833	work[kRow]=-1.0;
1834	}
1835	CoinBigIndex j;
1836	int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1837	for (j=0;j<number;j++) {
1838	int jRow = which[j];
1839	put[j]=work[jRow];
1840	work[jRow]=0.0;
1841	}
1842	}
1843	} else {
1844	// quadratic
1845	const int * columnQuadratic = quadratic->getIndices();
1846	const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1847	const int * columnQuadraticLength = quadratic->getVectorLengths();
1848	const double * quadraticElement = quadratic->getElements();
1849	for (iRow=0;iRow<numberRows_;iRow++) {
1850	longDouble * put = sparseFactor_+choleskyStart_[iRow];
1851	int * which = choleskyRow_+indexStart_[iRow];
1852	int iOriginalRow = permute_[iRow];
1853	if (iOriginalRow<numberColumns) {
1854	CoinBigIndex j;
1855	iColumn=iOriginalRow;
1856	longDouble value = diagonal[iColumn];
1857	if (fabs(value)>1.0e-100) {
1858	value = 1.0/value;
1859	for (j=columnQuadraticStart[iColumn];
1860	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1861	int jColumn = columnQuadratic[j];
1862	int jNewColumn = permuteInverse_[jColumn];
1863	if (jNewColumn>iRow) {
1864	work[jNewColumn]=-quadraticElement[j];
1865	} else if (iColumn==jColumn) {
1866	value += quadraticElement[j];
1867	}
1868	}
1869	largest = CoinMax(largest,fabs(value));
1870	diagonal_[iRow] = -value;
1871	CoinBigIndex start=columnStart[iColumn];
1872	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1873	for (j=start;j<end;j++) {
1874	int kRow = row[j]+numberTotal;
1875	kRow=permuteInverse_[kRow];
1876	if (kRow>iRow) {
1877	work[kRow]=element[j];
1878	largest = CoinMax(largest,fabs(element[j]));
1879	}
1880	}
1881	} else {
1882	diagonal_[iRow] = -value;
1883	}
1884	} else if (iOriginalRow<numberTotal) {
1885	longDouble value = diagonal[iOriginalRow];
1886	if (fabs(value)>1.0e-100) {
1887	value = 1.0/value;
1888	largest = CoinMax(largest,fabs(value));
1889	} else {
1890	value = 1.0e100;
1891	}
1892	diagonal_[iRow] = -value;
1893	int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
1894	if (kRow>iRow)
1895	work[kRow]=-1.0;
1896	} else {
1897	diagonal_[iRow]=delta2;
1898	int kRow = iOriginalRow-numberTotal;
1899	CoinBigIndex start=rowStart[kRow];
1900	CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
1901	for (CoinBigIndex j=start;j<end;j++) {
1902	int jRow = column[j];
1903	int jNewRow = permuteInverse_[jRow];
1904	if (jNewRow>iRow) {
1905	work[jNewRow]=elementByRow[j];
1906	largest = CoinMax(largest,fabs(elementByRow[j]));
1907	}
1908	}
1909	// slack - should it be permute
1910	kRow = permuteInverse_[kRow+numberColumns];
1911	if (kRow>iRow)
1912	work[kRow]=-1.0;
1913	}
1914	CoinBigIndex j;
1915	int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1916	for (j=0;j<number;j++) {
1917	int jRow = which[j];
1918	put[j]=work[jRow];
1919	work[jRow]=0.0;
1920	}
1921	for (j=0;j<numberRows_;j++)
1922	assert (!work[j]);
1923	}
1924	}
1925	} else {
1926	if (!quadratic) {
1927	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1928	longDouble * put = sparseFactor_+choleskyStart_[iColumn];
1929	int * which = choleskyRow_+indexStart_[iColumn];
1930	longDouble value = diagonal[iColumn];
1931	if (fabs(value)>1.0e-100) {
1932	value = 1.0/value;
1933	largest = CoinMax(largest,fabs(value));
1934	diagonal_[iColumn] = -value;
1935	CoinBigIndex start=columnStart[iColumn];
1936	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1937	for (CoinBigIndex j=start;j<end;j++) {
1938	//choleskyRow_[numberElements]=row[j]+numberTotal;
1939	//sparseFactor_[numberElements++]=element[j];
1940	work[row[j]+numberTotal]=element[j];
1941	largest = CoinMax(largest,fabs(element[j]));
1942	}
1943	} else {
1944	diagonal_[iColumn] = -value;
1945	}
1946	CoinBigIndex j;
1947	int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
1948	for (j=0;j<number;j++) {
1949	int jRow = which[j];
1950	put[j]=work[jRow];
1951	work[jRow]=0.0;
1952	}
1953	}
1954	} else {
1955	// Quadratic
1956	const int * columnQuadratic = quadratic->getIndices();
1957	const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1958	const int * columnQuadraticLength = quadratic->getVectorLengths();
1959	const double * quadraticElement = quadratic->getElements();
1960	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1961	longDouble * put = sparseFactor_+choleskyStart_[iColumn];
1962	int * which = choleskyRow_+indexStart_[iColumn];
1963	int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
1964	longDouble value = diagonal[iColumn];
1965	CoinBigIndex j;
1966	if (fabs(value)>1.0e-100) {
1967	value = 1.0/value;
1968	for (j=columnQuadraticStart[iColumn];
1969	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1970	int jColumn = columnQuadratic[j];
1971	if (jColumn>iColumn) {
1972	work[jColumn]=-quadraticElement[j];
1973	} else if (iColumn==jColumn) {
1974	value += quadraticElement[j];
1975	}
1976	}
1977	largest = CoinMax(largest,fabs(value));
1978	diagonal_[iColumn] = -value;
1979	CoinBigIndex start=columnStart[iColumn];
1980	CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1981	for (j=start;j<end;j++) {
1982	work[row[j]+numberTotal]=element[j];
1983	largest = CoinMax(largest,fabs(element[j]));
1984	}
1985	for (j=0;j<number;j++) {
1986	int jRow = which[j];
1987	put[j]=work[jRow];
1988	work[jRow]=0.0;
1989	}
1990	} else {
1991	value = 1.0e100;
1992	diagonal_[iColumn] = -value;
1993	for (j=0;j<number;j++) {
1994	int jRow = which[j];
1995	put[j]=work[jRow];
1996	}
1997	}
1998	}
1999	}
2000	// slacks
2001	for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
2002	longDouble * put = sparseFactor_+choleskyStart_[iColumn];
2003	int * which = choleskyRow_+indexStart_[iColumn];
2004	longDouble value = diagonal[iColumn];
2005	if (fabs(value)>1.0e-100) {
2006	value = 1.0/value;
2007	largest = CoinMax(largest,fabs(value));
2008	} else {
2009	value = 1.0e100;
2010	}
2011	diagonal_[iColumn] = -value;
2012	work[iColumn-numberColumns+numberTotal]=-1.0;
2013	CoinBigIndex j;
2014	int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
2015	for (j=0;j<number;j++) {
2016	int jRow = which[j];
2017	put[j]=work[jRow];
2018	work[jRow]=0.0;
2019	}
2020	}
2021	// Finish diagonal
2022	longDouble delta2 = model_->delta(); // add delta*delta to bottom
2023	delta2 *= delta2;
2024	for (iRow=0;iRow<numberRowsModel;iRow++) {
2025	longDouble * put = sparseFactor_+choleskyStart_[iRow+numberTotal];
2026	diagonal_[iRow+numberTotal]=delta2;
2027	CoinBigIndex j;
2028	int number = choleskyStart_[iRow+numberTotal+1]-choleskyStart_[iRow+numberTotal];
2029	for (j=0;j<number;j++) {
2030	put[j]=0.0;
2031	}
2032	}
2033	}
2034	//check sizes
2035	largest*=1.0e-20;
2036	largest = CoinMin(largest,1.0e-11);
2037	doubleParameters_[10]=CoinMax(1.0e-20,largest);
2038	integerParameters_[20]=0;
2039	doubleParameters_[3]=0.0;
2040	doubleParameters_[4]=COIN_DBL_MAX;
2041	// Set up LDL cutoff
2042	integerParameters_[34]=numberTotal;
2043	// KKT
2044	int * rowsDropped2 = new int[numberRows_];
2045	CoinZeroN(rowsDropped2,numberRows_);
2046	factorizePart2(rowsDropped2);
2047	newDropped=integerParameters_[20];
2048	largest=doubleParameters_[3];
2049	smallest=doubleParameters_[4];
2050	if (model_->messageHandler()->logLevel()>1)
2051	std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
2052	choleskyCondition_=largest/smallest;
2053	// Should save adjustments in ..R_
2054	int n1=0,n2=0;
2055	double * primalR = model_->primalR();
2056	double * dualR = model_->dualR();
2057	for (iRow=0;iRow<numberTotal;iRow++) {
2058	if (rowsDropped2[iRow]) {
2059	n1++;
2060	//printf("row region1 %d dropped\n",iRow);
2061	//rowsDropped_[iRow]=1;
2062	rowsDropped_[iRow]=0;
2063	primalR[iRow]=doubleParameters_[20];
2064	} else {
2065	rowsDropped_[iRow]=0;
2066	primalR[iRow]=0.0;
2067	}
2068	}
2069	for (;iRow<numberRows_;iRow++) {
2070	if (rowsDropped2[iRow]) {
2071	n2++;
2072	//printf("row region2 %d dropped\n",iRow);
2073	//rowsDropped_[iRow]=1;
2074	rowsDropped_[iRow]=0;
2075	dualR[iRow-numberTotal]=doubleParameters_[34];
2076	} else {
2077	rowsDropped_[iRow]=0;
2078	dualR[iRow-numberTotal]=0.0;
2079	}
2080	}
2081	delete [] rowsDropped2;
2082	}
2083	status_=0;
2084	return newDropped;
2085	}
2086	/* Factorize - filling in rowsDropped and returning number dropped
2087	in integerParam.
2088	*/
2089	void
2090	ClpCholeskyBase::factorizePart2(int * rowsDropped)
2091	{
2092	longDouble largest=doubleParameters_[3];
2093	longDouble smallest=doubleParameters_[4];
2094	int numberDropped=integerParameters_[20];
2095	// probably done before
2096	largest=0.0;
2097	smallest=COIN_DBL_MAX;
2098	numberDropped=0;
2099	double dropValue = doubleParameters_[10];
2100	int firstPositive=integerParameters_[34];
2101	longDouble * d = ClpCopyOfArray(diagonal_,numberRows_);
2102	int iRow;
2103	// minimum size before clique done
2104	//#define MINCLIQUE INT_MAX
2105	#define MINCLIQUE 3
2106	longDouble * work = workDouble_;
2107	CoinBigIndex * first = workInteger_;
2108
2109	for (iRow=0;iRow<numberRows_;iRow++) {
2110	link_[iRow]=-1;
2111	work[iRow]=0.0;
2112	first[iRow]=choleskyStart_[iRow];
2113	}
2114
2115	int lastClique=-1;
2116	bool inClique=false;
2117	bool newClique=false;
2118	bool endClique=false;
2119	int lastRow=0;
2120	int cliqueSize=0;
2121	CoinBigIndex cliquePointer=0;
2122	int nextRow2=-1;
2123
2124	for (iRow=0;iRow<firstDense_+1;iRow++) {
2125	if (iRow<firstDense_) {
2126	endClique=false;
2127	if (clique_[iRow]>0) {
2128	// this is in a clique
2129	inClique=true;
2130	if (clique_[iRow]>lastClique) {
2131	// new Clique
2132	newClique=true;
2133	// If we have clique going then signal to do old one
2134	endClique=(lastClique>0);
2135	} else {
2136	// Still in clique
2137	newClique=false;
2138	}
2139	} else {
2140	// not in clique
2141	inClique=false;
2142	newClique=false;
2143	// If we have clique going then signal to do old one
2144	endClique=(lastClique>0);
2145	}
2146	lastClique=clique_[iRow];
2147	} else if (inClique) {
2148	// Finish off
2149	endClique=true;
2150	} else {
2151	break;
2152	}
2153	if (endClique) {
2154	// We have just finished updating a clique - do block pivot and clean up
2155	int jRow;
2156	for ( jRow=lastRow;jRow<iRow;jRow++) {
2157	int jCount = jRow-lastRow;
2158	longDouble diagonalValue = diagonal_[jRow];
2159	CoinBigIndex start=choleskyStart_[jRow];
2160	CoinBigIndex end=choleskyStart_[jRow+1];
2161	for (int kRow=lastRow;kRow<jRow;kRow++) {
2162	jCount--;
2163	CoinBigIndex get = choleskyStart_[kRow]+jCount;
2164	longDouble a_jk = sparseFactor_[get];
2165	longDouble value1 = d[kRow]*a_jk;
2166	diagonalValue -= a_jk*value1;
2167	for (CoinBigIndex j=start;j<end;j++)
2168	sparseFactor_[j] -= value1*sparseFactor_[++get];
2169	}
2170	// check
2171	int originalRow = permute_[jRow];
2172	if (originalRow<firstPositive) {
2173	// must be negative
2174	if (diagonalValue<=-dropValue) {
2175	smallest = CoinMin(smallest,-diagonalValue);
2176	largest = CoinMax(largest,-diagonalValue);
2177	d[jRow]=diagonalValue;
2178	diagonalValue = 1.0/diagonalValue;
2179	} else {
2180	rowsDropped[originalRow]=2;
2181	d[jRow]=-1.0e100;
2182	diagonalValue=0.0;
2183	integerParameters_[20]++;
2184	}
2185	} else {
2186	// must be positive
2187	if (diagonalValue>=dropValue) {
2188	smallest = CoinMin(smallest,diagonalValue);
2189	largest = CoinMax(largest,diagonalValue);
2190	d[jRow]=diagonalValue;
2191	diagonalValue = 1.0/diagonalValue;
2192	} else {
2193	rowsDropped[originalRow]=2;
2194	d[jRow]=1.0e100;
2195	diagonalValue=0.0;
2196	integerParameters_[20]++;
2197	}
2198	}
2199	diagonal_[jRow]=diagonalValue;
2200	for (CoinBigIndex j=start;j<end;j++) {
2201	sparseFactor_[j] *= diagonalValue;
2202	}
2203	}
2204	if (nextRow2>=0) {
2205	for ( jRow=lastRow;jRow<iRow-1;jRow++) {
2206	link_[jRow]=jRow+1;
2207	}
2208	link_[iRow-1]=link_[nextRow2];
2209	link_[nextRow2]=lastRow;
2210	}
2211	}
2212	if (iRow==firstDense_)
2213	break; // we were just cleaning up
2214	if (newClique) {
2215	// initialize new clique
2216	lastRow=iRow;
2217	cliquePointer=choleskyStart_[iRow];
2218	cliqueSize=choleskyStart_[iRow+1]-cliquePointer+1;
2219	}
2220	// for each column L[,kRow] that affects L[,iRow]
2221	longDouble diagonalValue=diagonal_[iRow];
2222	int nextRow = link_[iRow];
2223	int kRow=0;
2224	while (1) {
2225	kRow=nextRow;
2226	if (kRow<0)
2227	break; // out of loop
2228	nextRow=link_[kRow];
2229	// Modify by outer product of L[,irow] by L[,krow] from first
2230	CoinBigIndex k=first[kRow];
2231	CoinBigIndex end=choleskyStart_[kRow+1];
2232	assert(k<end);
2233	longDouble a_ik=sparseFactor_[k++];
2234	longDouble value1=d[kRow]*a_ik;
2235	// update first
2236	first[kRow]=k;
2237	diagonalValue -= value1*a_ik;
2238	CoinBigIndex offset = indexStart_[kRow]-choleskyStart_[kRow];
2239	if (k<end) {
2240	int jRow = choleskyRow_[k+offset];
2241	if (clique_[kRow]<MINCLIQUE) {
2242	link_[kRow]=link_[jRow];
2243	link_[jRow]=kRow;
2244	for (;k<end;k++) {
2245	int jRow = choleskyRow_[k+offset];
2246	work[jRow] += sparseFactor_[k]*value1;
2247	}
2248	} else {
2249	// Clique
2250	CoinBigIndex currentIndex = k+offset;
2251	int linkSave=link_[jRow];
2252	link_[jRow]=kRow;
2253	work[kRow]=value1; // ? or a_jk
2254	int last = kRow+clique_[kRow];
2255	for (int kkRow=kRow+1;kkRow<last;kkRow++) {
2256	CoinBigIndex j=first[kkRow];
2257	//int iiRow = choleskyRow_[j+indexStart_[kkRow]-choleskyStart_[kkRow]];
2258	longDouble a = sparseFactor_[j];
2259	longDouble dValue = d[kkRow]*a;
2260	diagonalValue -= a*dValue;
2261	work[kkRow]=dValue;
2262	first[kkRow]++;
2263	link_[kkRow-1]=kkRow;
2264	}
2265	nextRow = link_[last-1];
2266	link_[last-1]=linkSave;
2267	int length = end-k;
2268	for (int i=0;i<length;i++) {
2269	int lRow = choleskyRow_[currentIndex++];
2270	longDouble t0 = work[lRow];
2271	for (int kkRow=kRow;kkRow<last;kkRow++) {
2272	CoinBigIndex j = first[kkRow]+i;
2273	t0 += work[kkRow]*sparseFactor_[j];
2274	}
2275	work[lRow]=t0;
2276	}
2277	}
2278	}
2279	}
2280	// Now apply
2281	if (inClique) {
2282	// in clique
2283	diagonal_[iRow]=diagonalValue;
2284	CoinBigIndex start=choleskyStart_[iRow];
2285	CoinBigIndex end=choleskyStart_[iRow+1];
2286	CoinBigIndex currentIndex = indexStart_[iRow];
2287	nextRow2=-1;
2288	CoinBigIndex get=start+clique_[iRow]-1;
2289	if (get<end) {
2290	nextRow2 = choleskyRow_[currentIndex+get-start];
2291	first[iRow]=get;
2292	}
2293	for (CoinBigIndex j=start;j<end;j++) {
2294	int kRow = choleskyRow_[currentIndex++];
2295	sparseFactor_[j] -= work[kRow]; // times?
2296	work[kRow]=0.0;
2297	}
2298	} else {
2299	// not in clique
2300	int originalRow = permute_[iRow];
2301	if (originalRow<firstPositive) {
2302	// must be negative
2303	if (diagonalValue<=-dropValue) {
2304	smallest = CoinMin(smallest,-diagonalValue);
2305	largest = CoinMax(largest,-diagonalValue);
2306	d[iRow]=diagonalValue;
2307	diagonalValue = 1.0/diagonalValue;
2308	} else {
2309	rowsDropped[originalRow]=2;
2310	d[iRow]=-1.0e100;
2311	diagonalValue=0.0;
2312	integerParameters_[20]++;
2313	}
2314	} else {
2315	// must be positive
2316	if (diagonalValue>=dropValue) {
2317	smallest = CoinMin(smallest,diagonalValue);
2318	largest = CoinMax(largest,diagonalValue);
2319	d[iRow]=diagonalValue;
2320	diagonalValue = 1.0/diagonalValue;
2321	} else {
2322	rowsDropped[originalRow]=2;
2323	d[iRow]=1.0e100;
2324	diagonalValue=0.0;
2325	integerParameters_[20]++;
2326	}
2327	}
2328	diagonal_[iRow]=diagonalValue;
2329	CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
2330	CoinBigIndex start = choleskyStart_[iRow];
2331	CoinBigIndex end = choleskyStart_[iRow+1];
2332	assert (first[iRow]==start);
2333	if (start<end) {
2334	int nextRow = choleskyRow_[start+offset];
2335	link_[iRow]=link_[nextRow];
2336	link_[nextRow]=iRow;
2337	for (CoinBigIndex j=start;j<end;j++) {
2338	int jRow = choleskyRow_[j+offset];
2339	longDouble value = sparseFactor_[j]-work[jRow];
2340	work[jRow]=0.0;
2341	sparseFactor_[j]= diagonalValue*value;
2342	}
2343	}
2344	}
2345	}
2346	if (firstDense_<numberRows_) {
2347	// do dense
2348	// update dense part
2349	updateDense(d,work,first);
2350	ClpCholeskyDense dense;
2351	// just borrow space
2352	int nDense = numberRows_-firstDense_;
2353	if (doKKT_) {
2354	for (iRow=firstDense_;iRow<numberRows_;iRow++) {
2355	int originalRow=permute_[iRow];
2356	if (originalRow>=firstPositive) {
2357	firstPositive=iRow-firstDense_;
2358	break;
2359	}
2360	}
2361	}
2362	dense.reserveSpace(this,nDense);
2363	int * dropped = new int[nDense];
2364	memset(dropped,0,nDense*sizeof(int));
2365	dense.setDoubleParameter(3,largest);
2366	dense.setDoubleParameter(4,smallest);
2367	dense.setDoubleParameter(10,dropValue);
2368	dense.setIntegerParameter(20,0);
2369	dense.setIntegerParameter(34,firstPositive);
2370	dense.factorizePart2(dropped);
2371	largest=dense.getDoubleParameter(3);
2372	smallest=dense.getDoubleParameter(4);
2373	integerParameters_[20]+=dense.getIntegerParameter(20);
2374	for (iRow=firstDense_;iRow<numberRows_;iRow++) {
2375	int originalRow=permute_[iRow];
2376	rowsDropped[originalRow]=dropped[iRow-firstDense_];
2377	}
2378	delete [] dropped;
2379	}
2380	delete [] d;
2381	doubleParameters_[3]=largest;
2382	doubleParameters_[4]=smallest;
2383	return;
2384	}
2385	// Updates dense part (broken out for profiling)
2386	void ClpCholeskyBase::updateDense(longDouble * d, longDouble * work, int * first)
2387	{
2388	for (int iRow=0;iRow<firstDense_;iRow++) {
2389	CoinBigIndex start=first[iRow];
2390	CoinBigIndex end=choleskyStart_[iRow+1];
2391	if (start<end) {
2392	CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
2393	if (clique_[iRow]<2) {
2394	longDouble dValue=d[iRow];
2395	for (CoinBigIndex k=start;k<end;k++) {
2396	int kRow = choleskyRow_[k+offset];
2397	assert(kRow>=firstDense_);
2398	longDouble a_ik=sparseFactor_[k];
2399	longDouble value1=dValue*a_ik;
2400	diagonal_[kRow] -= value1*a_ik;
2401	CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2402	for (CoinBigIndex j=k+1;j<end;j++) {
2403	int jRow=choleskyRow_[j+offset];
2404	longDouble a_jk = sparseFactor_[j];
2405	sparseFactor_[base+jRow] -= a_jk*value1;
2406	}
2407	}
2408	} else if (clique_[iRow]<3) {
2409	// do as pair
2410	longDouble dValue0=d[iRow];
2411	longDouble dValue1=d[iRow+1];
2412	int offset1 = first[iRow+1]-start;
2413	// skip row
2414	iRow++;
2415	for (CoinBigIndex k=start;k<end;k++) {
2416	int kRow = choleskyRow_[k+offset];
2417	assert(kRow>=firstDense_);
2418	longDouble a_ik0=sparseFactor_[k];
2419	longDouble value0=dValue0*a_ik0;
2420	longDouble a_ik1=sparseFactor_[k+offset1];
2421	longDouble value1=dValue1*a_ik1;
2422	diagonal_[kRow] -= value0a_ik0 + value1a_ik1;
2423	CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2424	for (CoinBigIndex j=k+1;j<end;j++) {
2425	int jRow=choleskyRow_[j+offset];
2426	longDouble a_jk0 = sparseFactor_[j];
2427	longDouble a_jk1 = sparseFactor_[j+offset1];
2428	sparseFactor_[base+jRow] -= a_jk0value0+a_jk1value1;
2429	}
2430	}
2431	#define MANY_REGISTERS
2432	#ifdef MANY_REGISTERS
2433	} else if (clique_[iRow]==3) {
2434	#else
2435	} else {
2436	#endif
2437	// do as clique
2438	// maybe later get fancy on big cliques and do transpose copy
2439	// seems only worth going to 3 on Intel
2440	longDouble dValue0=d[iRow];
2441	longDouble dValue1=d[iRow+1];
2442	longDouble dValue2=d[iRow+2];
2443	// get offsets and skip rows
2444	int offset1 = first[++iRow]-start;
2445	int offset2 = first[++iRow]-start;
2446	for (CoinBigIndex k=start;k<end;k++) {
2447	int kRow = choleskyRow_[k+offset];
2448	assert(kRow>=firstDense_);
2449	double diagonalValue=diagonal_[kRow];
2450	longDouble a_ik0=sparseFactor_[k];
2451	longDouble value0=dValue0*a_ik0;
2452	longDouble a_ik1=sparseFactor_[k+offset1];
2453	longDouble value1=dValue1*a_ik1;
2454	longDouble a_ik2=sparseFactor_[k+offset2];
2455	longDouble value2=dValue2*a_ik2;
2456	CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2457	diagonal_[kRow] = diagonalValue - value0a_ik0 - value1a_ik1 - value2*a_ik2;
2458	for (CoinBigIndex j=k+1;j<end;j++) {
2459	int jRow=choleskyRow_[j+offset];
2460	longDouble a_jk0 = sparseFactor_[j];
2461	longDouble a_jk1 = sparseFactor_[j+offset1];
2462	longDouble a_jk2 = sparseFactor_[j+offset2];
2463	sparseFactor_[base+jRow] -= a_jk0value0+a_jk1value1+a_jk2*value2;
2464	}
2465	}
2466	#ifdef MANY_REGISTERS
2467	} else {
2468	// do as clique
2469	// maybe later get fancy on big cliques and do transpose copy
2470	// maybe only worth going to 3 on Intel (but may have hidden registers)
2471	longDouble dValue0=d[iRow];
2472	longDouble dValue1=d[iRow+1];
2473	longDouble dValue2=d[iRow+2];
2474	longDouble dValue3=d[iRow+3];
2475	// get offsets and skip rows
2476	int offset1 = first[++iRow]-start;
2477	int offset2 = first[++iRow]-start;
2478	int offset3 = first[++iRow]-start;
2479	for (CoinBigIndex k=start;k<end;k++) {
2480	int kRow = choleskyRow_[k+offset];
2481	assert(kRow>=firstDense_);
2482	double diagonalValue=diagonal_[kRow];
2483	longDouble a_ik0=sparseFactor_[k];
2484	longDouble value0=dValue0*a_ik0;
2485	longDouble a_ik1=sparseFactor_[k+offset1];
2486	longDouble value1=dValue1*a_ik1;
2487	longDouble a_ik2=sparseFactor_[k+offset2];
2488	longDouble value2=dValue2*a_ik2;
2489	longDouble a_ik3=sparseFactor_[k+offset3];
2490	longDouble value3=dValue3*a_ik3;
2491	CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2492	diagonal_[kRow] = diagonalValue - (value0a_ik0 + value1a_ik1 + value2a_ik2+value3a_ik3);
2493	for (CoinBigIndex j=k+1;j<end;j++) {
2494	int jRow=choleskyRow_[j+offset];
2495	longDouble a_jk0 = sparseFactor_[j];
2496	longDouble a_jk1 = sparseFactor_[j+offset1];
2497	longDouble a_jk2 = sparseFactor_[j+offset2];
2498	longDouble a_jk3 = sparseFactor_[j+offset3];
2499	sparseFactor_[base+jRow] -= a_jk0value0+a_jk1value1+a_jk2value2+a_jk3value3;
2500	}
2501	}
2502	#endif
2503	}
2504	}
2505	}
2506	}
2507	/* Uses factorization to solve. */
2508	void
2509	ClpCholeskyBase::solve (double * region)
2510	{
2511	if (!whichDense_) {
2512	solve(region,3);
2513	} else {
2514	// dense columns
2515	int i;
2516	solve(region,1);
2517	// do change;
2518	int numberDense = dense_->numberRows();
2519	double * change = new double[numberDense];
2520	for (i=0;i<numberDense;i++) {
2521	const longDouble * a = denseColumn_+i*numberRows_;
2522	longDouble value =0.0;
2523	for (int iRow=0;iRow<numberRows_;iRow++)
2524	value += a[iRow]*region[iRow];
2525	change[i]=value;
2526	}
2527	// solve
2528	dense_->solve(change);
2529	for (i=0;i<numberDense;i++) {
2530	const longDouble * a = denseColumn_+i*numberRows_;
2531	longDouble value = change[i];
2532	for (int iRow=0;iRow<numberRows_;iRow++)
2533	region[iRow] -= value*a[iRow];
2534	}
2535	delete [] change;
2536	// and finish off
2537	solve(region,2);
2538	}
2539	}
2540	/* solve - 1 just first half, 2 just second half - 3 both.
2541	If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
2542	*/
2543	void
2544	ClpCholeskyBase::solve(double * region, int type)
2545	{
2546	#ifdef CLP_DEBUG
2547	double * regionX=NULL;
2548	if (sizeof(longWork)!=sizeof(double)&&type==3) {
2549	regionX=ClpCopyOfArray(region,numberRows_);
2550	}
2551	#endif
2552	longWork * work = (longWork *) workDouble_;
2553	int i;
2554	CoinBigIndex j;
2555	for (i=0;i<numberRows_;i++) {
2556	int iRow = permute_[i];
2557	work[i] = region[iRow];
2558	}
2559	switch (type) {
2560	case 1:
2561	for (i=0;i<numberRows_;i++) {
2562	longDouble value=work[i];
2563	CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2564	for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2565	int iRow = choleskyRow_[j+offset];
2566	work[iRow] -= sparseFactor_[j]*value;
2567	}
2568	}
2569	for (i=0;i<numberRows_;i++) {
2570	int iRow = permute_[i];
2571	region[iRow]=work[i]*diagonal_[i];
2572	}
2573	break;
2574	case 2:
2575	for (i=numberRows_-1;i>=0;i--) {
2576	CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2577	longDouble value=work[i]*diagonal_[i];
2578	for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2579	int iRow = choleskyRow_[j+offset];
2580	value -= sparseFactor_[j]*work[iRow];
2581	}
2582	work[i]=value;
2583	int iRow = permute_[i];
2584	region[iRow]=value;
2585	}
2586	break;
2587	case 3:
2588	for (i=0;i<firstDense_;i++) {
2589	CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2590	longDouble value=work[i];
2591	for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2592	int iRow = choleskyRow_[j+offset];
2593	work[iRow] -= sparseFactor_[j]*value;
2594	}
2595	}
2596	if (firstDense_<numberRows_) {
2597	// do dense
2598	ClpCholeskyDense dense;
2599	// just borrow space
2600	int nDense = numberRows_-firstDense_;
2601	dense.reserveSpace(this,nDense);
2602	#if CLP_LONG_CHOLESKY!=1
2603	dense.solveLong(work+firstDense_);
2604	#else
2605	dense.solveLongWork(work+firstDense_);
2606	#endif
2607	for (i=numberRows_-1;i>=firstDense_;i--) {
2608	longDouble value=work[i];
2609	int iRow = permute_[i];
2610	region[iRow]=value;
2611	}
2612	}
2613	for (i=firstDense_-1;i>=0;i--) {
2614	CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2615	longDouble value=work[i]*diagonal_[i];
2616	for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2617	int iRow = choleskyRow_[j+offset];
2618	value -= sparseFactor_[j]*work[iRow];
2619	}
2620	work[i]=value;
2621	int iRow = permute_[i];
2622	region[iRow]=value;
2623	}
2624	break;
2625	}
2626	#ifdef CLP_DEBUG
2627	if (regionX) {
2628	longDouble * work = workDouble_;
2629	int i;
2630	CoinBigIndex j;
2631	double largestO=0.0;
2632	for (i=0;i<numberRows_;i++) {
2633	largestO = CoinMax(largestO,fabs(regionX[i]));
2634	}
2635	for (i=0;i<numberRows_;i++) {
2636	int iRow = permute_[i];
2637	work[i] = regionX[iRow];
2638	}
2639	for (i=0;i<firstDense_;i++) {
2640	CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2641	longDouble value=work[i];
2642	for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2643	int iRow = choleskyRow_[j+offset];
2644	work[iRow] -= sparseFactor_[j]*value;
2645	}
2646	}
2647	if (firstDense_<numberRows_) {
2648	// do dense
2649	ClpCholeskyDense dense;
2650	// just borrow space
2651	int nDense = numberRows_-firstDense_;
2652	dense.reserveSpace(this,nDense);
2653	dense.solveLong(work+firstDense_);
2654	for (i=numberRows_-1;i>=firstDense_;i--) {
2655	longDouble value=work[i];
2656	int iRow = permute_[i];
2657	regionX[iRow]=value;
2658	}
2659	}
2660	for (i=firstDense_-1;i>=0;i--) {
2661	CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2662	longDouble value=work[i]*diagonal_[i];
2663	for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2664	int iRow = choleskyRow_[j+offset];
2665	value -= sparseFactor_[j]*work[iRow];
2666	}
2667	work[i]=value;
2668	int iRow = permute_[i];
2669	regionX[iRow]=value;
2670	}
2671	double largest=0.0;
2672	double largestV=0.0;
2673	for (i=0;i<numberRows_;i++) {
2674	largest = CoinMax(largest,fabs(region[i]-regionX[i]));
2675	largestV = CoinMax(largestV,fabs(region[i]));
2676	}
2677	printf("largest difference %g, largest %g, largest original %g\n",
2678	largest,largestV,largestO);
2679	delete [] regionX;
2680	}
2681	#endif
2682	}
2683	/* solve - 1 just first half, 2 just second half - 3 both.
2684	If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
2685	*/
2686	void
2687	ClpCholeskyBase::solveLong(longDouble * region, int type)
2688	{
2689	int i;
2690	CoinBigIndex j;
2691	for (i=0;i<numberRows_;i++) {
2692	int iRow = permute_[i];
2693	workDouble_[i] = region[iRow];
2694	}
2695	switch (type) {
2696	case 1:
2697	for (i=0;i<numberRows_;i++) {
2698	longDouble value=workDouble_[i];
2699	CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2700	for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2701	int iRow = choleskyRow_[j+offset];
2702	workDouble_[iRow] -= sparseFactor_[j]*value;
2703	}
2704	}
2705	for (i=0;i<numberRows_;i++) {
2706	int iRow = permute_[i];
2707	region[iRow]=workDouble_[i]*diagonal_[i];
2708	}
2709	break;
2710	case 2:
2711	for (i=numberRows_-1;i>=0;i--) {
2712	CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2713	longDouble value=workDouble_[i]*diagonal_[i];
2714	for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2715	int iRow = choleskyRow_[j+offset];
2716	value -= sparseFactor_[j]*workDouble_[iRow];
2717	}
2718	workDouble_[i]=value;
2719	int iRow = permute_[i];
2720	region[iRow]=value;
2721	}
2722	break;
2723	case 3:
2724	for (i=0;i<firstDense_;i++) {
2725	CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2726	longDouble value=workDouble_[i];
2727	for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2728	int iRow = choleskyRow_[j+offset];
2729	workDouble_[iRow] -= sparseFactor_[j]*value;
2730	}
2731	}
2732	if (firstDense_<numberRows_) {
2733	// do dense
2734	ClpCholeskyDense dense;
2735	// just borrow space
2736	int nDense = numberRows_-firstDense_;
2737	dense.reserveSpace(this,nDense);
2738	dense.solveLong(workDouble_+firstDense_);
2739	for (i=numberRows_-1;i>=firstDense_;i--) {
2740	longDouble value=workDouble_[i];
2741	int iRow = permute_[i];
2742	region[iRow]=value;
2743	}
2744	}
2745	for (i=firstDense_-1;i>=0;i--) {
2746	CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2747	longDouble value=workDouble_[i]*diagonal_[i];
2748	for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2749	int iRow = choleskyRow_[j+offset];
2750	value -= sparseFactor_[j]*workDouble_[iRow];
2751	}
2752	workDouble_[i]=value;
2753	int iRow = permute_[i];
2754	region[iRow]=value;
2755	}
2756	break;
2757	}
2758	}
2759

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

Context Navigation

source: trunk/ClpCholeskyBase.cpp @ 629

Download in other formats: