Context Navigation

source: trunk/Clp/src/ClpCholeskyMumps.cpp @ 1525

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago
Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.
Property svn:keywords set to `Id`
File size: 15.9 KB

Line
1	/* $Id: ClpCholeskyMumps.cpp 1525 2010-02-26 17:27:59Z mjs $ */
2	#ifdef MUMPS_BARRIER
3	// Copyright (C) 2004, International Business Machines
4	// Corporation and others. All Rights Reserved.
5
6
7
8	#include "CoinPragma.hpp"
9	#include "ClpCholeskyMumps.hpp"
10	extern "C"
11	{
12	//#include "mpi.h"
13	#include "/home/jjforre/cbc-trunk/ThirdParty/Mumps/MUMPS/libseq/mpi.h"
14	}
15	#define JOB_INIT -1
16	#define JOB_END -2
17	#define USE_COMM_WORLD -987654
18	#include "ClpMessage.hpp"
19	#include "ClpInterior.hpp"
20	#include "CoinHelperFunctions.hpp"
21	#include "ClpHelperFunctions.hpp"
22	//#############################################################################
23	// Constructors / Destructor / Assignment
24	//#############################################################################
25
26	//-------------------------------------------------------------------
27	// Default Constructor
28	//-------------------------------------------------------------------
29	ClpCholeskyMumps::ClpCholeskyMumps (int denseThreshold)
30	: ClpCholeskyBase(denseThreshold)
31	{
32	type_ = 16;
33	mumps_.n = 0;
34	mumps_.nz = 0;
35	mumps_.a = NULL;
36	mumps_.jcn = NULL;
37	mumps_.irn = NULL;
38	mumps_.job = JOB_INIT;//initialize mumps
39	mumps_.par = 1;//working host for sequential version
40	mumps_.sym = 2;//general symmetric matrix
41	mumps_.comm_fortran = USE_COMM_WORLD;
42	int myid, ierr;
43	int justName;
44	ierr = MPI_Init(&justName, NULL);
45	ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);
46	assert (!ierr);
47	dmumps_c(&mumps_);
48	#define ICNTL(I) icntl[(I)-1] /* macro s.t. indices match documentation */
49	#define CNTL(I) cntl[(I)-1] /* macro s.t. indices match documentation */
50	mumps_.ICNTL(5) = 1; // say compressed format
51	mumps_.ICNTL(4) = 2; // log messages
52	mumps_.ICNTL(24) = 1; // Deal with zeros on diagonal
53	mumps_.CNTL(3) = 1.0e-20; // drop if diagonal less than this
54	// output off
55	mumps_.ICNTL(1) = -1;
56	mumps_.ICNTL(2) = -1;
57	mumps_.ICNTL(3) = -1;
58	mumps_.ICNTL(4) = 0;
59	}
60
61	//-------------------------------------------------------------------
62	// Copy constructor
63	//-------------------------------------------------------------------
64	ClpCholeskyMumps::ClpCholeskyMumps (const ClpCholeskyMumps & rhs)
65	: ClpCholeskyBase(rhs)
66	{
67	abort();
68	}
69
70
71	//-------------------------------------------------------------------
72	// Destructor
73	//-------------------------------------------------------------------
74	ClpCholeskyMumps::~ClpCholeskyMumps ()
75	{
76	mumps_.job = JOB_END;
77	dmumps_c(&mumps_); /* Terminate instance */
78	MPI_Finalize();
79	}
80
81	//----------------------------------------------------------------
82	// Assignment operator
83	//-------------------------------------------------------------------
84	ClpCholeskyMumps &
85	ClpCholeskyMumps::operator=(const ClpCholeskyMumps& rhs)
86	{
87	if (this != &rhs) {
88	ClpCholeskyBase::operator=(rhs);
89	abort();
90	}
91	return *this;
92	}
93	//-------------------------------------------------------------------
94	// Clone
95	//-------------------------------------------------------------------
96	ClpCholeskyBase * ClpCholeskyMumps::clone() const
97	{
98	return new ClpCholeskyMumps(*this);
99	}
100	/* Orders rows and saves pointer to matrix.and model */
101	int
102	ClpCholeskyMumps::order(ClpInterior * model)
103	{
104	numberRows_ = model->numberRows();
105	if (doKKT_) {
106	numberRows_ += numberRows_ + model->numberColumns();
107	printf("finish coding MUMPS KKT!\n");
108	abort();
109	}
110	rowsDropped_ = new char [numberRows_];
111	memset(rowsDropped_, 0, numberRows_);
112	numberRowsDropped_ = 0;
113	model_ = model;
114	rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
115	const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
116	const int * columnLength = model_->clpMatrix()->getVectorLengths();
117	const int * row = model_->clpMatrix()->getIndices();
118	const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
119	const int * rowLength = rowCopy_->getVectorLengths();
120	const int * column = rowCopy_->getIndices();
121	// We need two arrays for counts
122	int * which = new int [numberRows_];
123	int * used = new int[numberRows_+1];
124	CoinZeroN(used, numberRows_);
125	int iRow;
126	sizeFactor_ = 0;
127	for (iRow = 0; iRow < numberRows_; iRow++) {
128	int number = 1;
129	// make sure diagonal exists
130	which[0] = iRow;
131	used[iRow] = 1;
132	if (!rowsDropped_[iRow]) {
133	CoinBigIndex startRow = rowStart[iRow];
134	CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
135	for (CoinBigIndex k = startRow; k < endRow; k++) {
136	int iColumn = column[k];
137	CoinBigIndex start = columnStart[iColumn];
138	CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
139	for (CoinBigIndex j = start; j < end; j++) {
140	int jRow = row[j];
141	if (jRow >= iRow && !rowsDropped_[jRow]) {
142	if (!used[jRow]) {
143	used[jRow] = 1;
144	which[number++] = jRow;
145	}
146	}
147	}
148	}
149	sizeFactor_ += number;
150	int j;
151	for (j = 0; j < number; j++)
152	used[which[j]] = 0;
153	}
154	}
155	delete [] which;
156	// NOT COMPRESSED FOR NOW ??? - Space for starts
157	mumps_.ICNTL(5) = 0; // say NOT compressed format
158	try {
159	choleskyStart_ = new CoinBigIndex[numberRows_+1+sizeFactor_];
160	} catch (...) {
161	// no memory
162	return -1;
163	}
164	// Now we have size - create arrays and fill in
165	try {
166	choleskyRow_ = new int [sizeFactor_];
167	} catch (...) {
168	// no memory
169	delete [] choleskyStart_;
170	choleskyStart_ = NULL;
171	return -1;
172	}
173	try {
174	sparseFactor_ = new double[sizeFactor_];
175	} catch (...) {
176	// no memory
177	delete [] choleskyRow_;
178	choleskyRow_ = NULL;
179	delete [] choleskyStart_;
180	choleskyStart_ = NULL;
181	return -1;
182	}
183
184	sizeFactor_ = 0;
185	which = choleskyRow_;
186	for (iRow = 0; iRow < numberRows_; iRow++) {
187	int number = 1;
188	// make sure diagonal exists
189	which[0] = iRow;
190	used[iRow] = 1;
191	choleskyStart_[iRow] = sizeFactor_;
192	if (!rowsDropped_[iRow]) {
193	CoinBigIndex startRow = rowStart[iRow];
194	CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
195	for (CoinBigIndex k = startRow; k < endRow; k++) {
196	int iColumn = column[k];
197	CoinBigIndex start = columnStart[iColumn];
198	CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
199	for (CoinBigIndex j = start; j < end; j++) {
200	int jRow = row[j];
201	if (jRow >= iRow && !rowsDropped_[jRow]) {
202	if (!used[jRow]) {
203	used[jRow] = 1;
204	which[number++] = jRow;
205	}
206	}
207	}
208	}
209	sizeFactor_ += number;
210	int j;
211	for (j = 0; j < number; j++)
212	used[which[j]] = 0;
213	// Sort
214	std::sort(which, which + number);
215	// move which on
216	which += number;
217	}
218	}
219	choleskyStart_[numberRows_] = sizeFactor_;
220	delete [] used;
221	permuteInverse_ = new int [numberRows_];
222	permute_ = new int[numberRows_];
223	// To Fortran and fake
224	for (iRow = 0; iRow < numberRows_ + 1; iRow++) {
225	int k = choleskyStart_[iRow];
226	int kEnd = choleskyStart_[iRow+1];
227	k += numberRows_ + 1;
228	kEnd += numberRows_ + 1;
229	for (; k < kEnd; k++)
230	choleskyStart_[k] = iRow + 1;
231	choleskyStart_[iRow]++;
232	}
233	mumps_.nz = sizeFactor_;
234	mumps_.irn = choleskyStart_ + numberRows_ + 1;
235	mumps_.jcn = choleskyRow_;
236	mumps_.a = NULL;
237	for (CoinBigIndex i = 0; i < sizeFactor_; i++) {
238	choleskyRow_[i]++;
239	#ifndef NDEBUG
240	assert (mumps_.irn[i] >= 1 && mumps_.irn[i] <= numberRows_);
241	assert (mumps_.jcn[i] >= 1 && mumps_.jcn[i] <= numberRows_);
242	#endif
243	}
244	// validate
245	//mumps code here
246	mumps_.n = numberRows_;
247	mumps_.nelt = numberRows_;
248	mumps_.eltptr = choleskyStart_;
249	mumps_.eltvar = choleskyRow_;
250	mumps_.a_elt = NULL;
251	mumps_.rhs = NULL;
252	mumps_.job = 1; // order
253	dmumps_c(&mumps_);
254	mumps_.a = sparseFactor_;
255	if (mumps_.infog[0]) {
256	printf("MUMPS ordering failed -error %d %d\n",
257	mumps_.infog[0], mumps_.infog[1]);
258	return 1;
259	} else {
260	double size = mumps_.infog[19];
261	if (size < 0)
262	size *= -1000000;
263	printf("%g nonzeros, flop count %g\n", size, mumps_.rinfog[0]);
264	}
265	for (iRow = 0; iRow < numberRows_; iRow++) {
266	permuteInverse_[iRow] = iRow;
267	permute_[iRow] = iRow;
268	}
269	return 0;
270	}
271	/* Does Symbolic factorization given permutation.
272	This is called immediately after order. If user provides this then
273	user must provide factorize and solve. Otherwise the default factorization is used
274	returns non-zero if not enough memory */
275	int
276	ClpCholeskyMumps::symbolic()
277	{
278	return 0;
279	}
280	/* Factorize - filling in rowsDropped and returning number dropped */
281	int
282	ClpCholeskyMumps::factorize(const double * diagonal, int * rowsDropped)
283	{
284	const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
285	const int * columnLength = model_->clpMatrix()->getVectorLengths();
286	const int * row = model_->clpMatrix()->getIndices();
287	const double * element = model_->clpMatrix()->getElements();
288	const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
289	const int * rowLength = rowCopy_->getVectorLengths();
290	const int * column = rowCopy_->getIndices();
291	const double * elementByRow = rowCopy_->getElements();
292	int numberColumns = model_->clpMatrix()->getNumCols();
293	int iRow;
294	double * work = new double[numberRows_];
295	CoinZeroN(work, numberRows_);
296	const double * diagonalSlack = diagonal + numberColumns;
297	int newDropped = 0;
298	double largest;
299	//double smallest;
300	//perturbation
301	double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
302	perturbation = 0.0;
303	perturbation = perturbation * perturbation;
304	if (perturbation > 1.0) {
305	#ifdef COIN_DEVELOP
306	//if (model_->model()->logLevel()&4)
307	std::cout << "large perturbation " << perturbation << std::endl;
308	#endif
309	perturbation = sqrt(perturbation);;
310	perturbation = 1.0;
311	}
312	double delta2 = model_->delta(); // add delta*delta to diagonal
313	delta2 *= delta2;
314	for (iRow = 0; iRow < numberRows_; iRow++) {
315	double * put = sparseFactor_ + choleskyStart_[iRow] - 1; // Fortran
316	int * which = choleskyRow_ + choleskyStart_[iRow] - 1; // Fortran
317	int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
318	if (!rowLength[iRow])
319	rowsDropped_[iRow] = 1;
320	if (!rowsDropped_[iRow]) {
321	CoinBigIndex startRow = rowStart[iRow];
322	CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
323	work[iRow] = diagonalSlack[iRow] + delta2;
324	for (CoinBigIndex k = startRow; k < endRow; k++) {
325	int iColumn = column[k];
326	if (!whichDense_ \|\| !whichDense_[iColumn]) {
327	CoinBigIndex start = columnStart[iColumn];
328	CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
329	double multiplier = diagonal[iColumn] * elementByRow[k];
330	for (CoinBigIndex j = start; j < end; j++) {
331	int jRow = row[j];
332	if (jRow >= iRow && !rowsDropped_[jRow]) {
333	double value = element[j] * multiplier;
334	work[jRow] += value;
335	}
336	}
337	}
338	}
339	int j;
340	for (j = 0; j < number; j++) {
341	int jRow = which[j] - 1; // to Fortran
342	put[j] = work[jRow];
343	work[jRow] = 0.0;
344	}
345	} else {
346	// dropped
347	int j;
348	for (j = 1; j < number; j++) {
349	put[j] = 0.0;
350	}
351	put[0] = 1.0;
352	}
353	}
354	//check sizes
355	double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_);
356	largest2 *= 1.0e-20;
357	largest = CoinMin(largest2, 1.0e-11);
358	int numberDroppedBefore = 0;
359	for (iRow = 0; iRow < numberRows_; iRow++) {
360	int dropped = rowsDropped_[iRow];
361	// Move to int array
362	rowsDropped[iRow] = dropped;
363	if (!dropped) {
364	CoinBigIndex start = choleskyStart_[iRow] - 1; // to Fortran
365	double diagonal = sparseFactor_[start];
366	if (diagonal > largest2) {
367	sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
368	} else {
369	sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
370	rowsDropped[iRow] = 2;
371	numberDroppedBefore++;
372	}
373	}
374	}
375	delete [] work;
376	// code here
377	mumps_.a_elt = sparseFactor_;
378	mumps_.rhs = NULL;
379	mumps_.job = 2; // factorize
380	dmumps_c(&mumps_);
381	if (mumps_.infog[0]) {
382	printf("MUMPS factorization failed -error %d %d\n",
383	mumps_.infog[0], mumps_.infog[1]);
384	}
385	choleskyCondition_ = 1.0;
386	bool cleanCholesky;
387	if (model_->numberIterations() < 2000)
388	cleanCholesky = true;
389	else
390	cleanCholesky = false;
391	if (cleanCholesky) {
392	//drop fresh makes some formADAT easier
393	//int oldDropped=numberRowsDropped_;
394	if (newDropped \|\| numberRowsDropped_) {
395	//std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
396	// newDropped<<" dropped)";
397	//if (newDropped>oldDropped)
398	//std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
399	//std::cout<<std::endl;
400	newDropped = 0;
401	for (int i = 0; i < numberRows_; i++) {
402	char dropped = rowsDropped[i];
403	rowsDropped_[i] = dropped;
404	if (dropped == 2) {
405	//dropped this time
406	rowsDropped[newDropped++] = i;
407	rowsDropped_[i] = 0;
408	}
409	}
410	numberRowsDropped_ = newDropped;
411	newDropped = -(2 + newDropped);
412	}
413	} else {
414	if (newDropped) {
415	newDropped = 0;
416	for (int i = 0; i < numberRows_; i++) {
417	char dropped = rowsDropped[i];
418	rowsDropped_[i] = dropped;
419	if (dropped == 2) {
420	//dropped this time
421	rowsDropped[newDropped++] = i;
422	rowsDropped_[i] = 1;
423	}
424	}
425	}
426	numberRowsDropped_ += newDropped;
427	if (numberRowsDropped_ && 0) {
428	std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " <<
429	numberRowsDropped_ << " dropped)";
430	if (newDropped) {
431	std::cout << " ( " << newDropped << " dropped this time)";
432	}
433	std::cout << std::endl;
434	}
435	}
436	status_ = 0;
437	return newDropped;
438	}
439	/* Uses factorization to solve. */
440	void
441	ClpCholeskyMumps::solve (double * region)
442	{
443	mumps_.rhs = region;
444	mumps_.job = 3; // solve
445	dmumps_c(&mumps_);
446	}
447	#endif

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

Download in other formats: