limp-cbc-0.3.2.0: cbits/coin/ClpHelperFunctions.cpp
/* $Id: ClpHelperFunctions.cpp 1817 2011-11-03 09:26:23Z forrest $ */
// Copyright (C) 2003, International Business Machines
// Corporation and others. All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).
/*
Note (JJF) I have added some operations on arrays even though they may
duplicate CoinDenseVector. I think the use of templates was a mistake
as I don't think inline generic code can take as much advantage of
parallelism or machine architectures or memory hierarchies.
*/
#include <cfloat>
#include <cstdlib>
#include <cmath>
#include "CoinHelperFunctions.hpp"
#include "CoinTypes.hpp"
double
maximumAbsElement(const double * region, int size)
{
int i;
double maxValue = 0.0;
for (i = 0; i < size; i++)
maxValue = CoinMax(maxValue, fabs(region[i]));
return maxValue;
}
void
setElements(double * region, int size, double value)
{
int i;
for (i = 0; i < size; i++)
region[i] = value;
}
void
multiplyAdd(const double * region1, int size, double multiplier1,
double * region2, double multiplier2)
{
int i;
if (multiplier1 == 1.0) {
if (multiplier2 == 1.0) {
for (i = 0; i < size; i++)
region2[i] = region1[i] + region2[i];
} else if (multiplier2 == -1.0) {
for (i = 0; i < size; i++)
region2[i] = region1[i] - region2[i];
} else if (multiplier2 == 0.0) {
for (i = 0; i < size; i++)
region2[i] = region1[i] ;
} else {
for (i = 0; i < size; i++)
region2[i] = region1[i] + multiplier2 * region2[i];
}
} else if (multiplier1 == -1.0) {
if (multiplier2 == 1.0) {
for (i = 0; i < size; i++)
region2[i] = -region1[i] + region2[i];
} else if (multiplier2 == -1.0) {
for (i = 0; i < size; i++)
region2[i] = -region1[i] - region2[i];
} else if (multiplier2 == 0.0) {
for (i = 0; i < size; i++)
region2[i] = -region1[i] ;
} else {
for (i = 0; i < size; i++)
region2[i] = -region1[i] + multiplier2 * region2[i];
}
} else if (multiplier1 == 0.0) {
if (multiplier2 == 1.0) {
// nothing to do
} else if (multiplier2 == -1.0) {
for (i = 0; i < size; i++)
region2[i] = -region2[i];
} else if (multiplier2 == 0.0) {
for (i = 0; i < size; i++)
region2[i] = 0.0;
} else {
for (i = 0; i < size; i++)
region2[i] = multiplier2 * region2[i];
}
} else {
if (multiplier2 == 1.0) {
for (i = 0; i < size; i++)
region2[i] = multiplier1 * region1[i] + region2[i];
} else if (multiplier2 == -1.0) {
for (i = 0; i < size; i++)
region2[i] = multiplier1 * region1[i] - region2[i];
} else if (multiplier2 == 0.0) {
for (i = 0; i < size; i++)
region2[i] = multiplier1 * region1[i] ;
} else {
for (i = 0; i < size; i++)
region2[i] = multiplier1 * region1[i] + multiplier2 * region2[i];
}
}
}
double
innerProduct(const double * region1, int size, const double * region2)
{
int i;
double value = 0.0;
for (i = 0; i < size; i++)
value += region1[i] * region2[i];
return value;
}
void
getNorms(const double * region, int size, double & norm1, double & norm2)
{
norm1 = 0.0;
norm2 = 0.0;
int i;
for (i = 0; i < size; i++) {
norm2 += region[i] * region[i];
norm1 = CoinMax(norm1, fabs(region[i]));
}
}
#ifndef NDEBUG
#include "ClpModel.hpp"
#include "ClpMessage.hpp"
ClpModel * clpTraceModel=NULL; // Set to trap messages
void ClpTracePrint(std::string fileName, std::string message, int lineNumber)
{
if (!clpTraceModel) {
std::cout<<fileName<<":"<<lineNumber<<" : \'"
<<message<<"\' failed."<<std::endl;
} else {
char line[1000];
sprintf(line,"%s: %d : \'%s\' failed.",fileName.c_str(),lineNumber,message.c_str());
clpTraceModel->messageHandler()->message(CLP_GENERAL_WARNING, clpTraceModel->messages())
<< line
<< CoinMessageEol;
}
}
#endif
#if COIN_LONG_WORK
// For long double versions
CoinWorkDouble
maximumAbsElement(const CoinWorkDouble * region, int size)
{
int i;
CoinWorkDouble maxValue = 0.0;
for (i = 0; i < size; i++)
maxValue = CoinMax(maxValue, CoinAbs(region[i]));
return maxValue;
}
void
setElements(CoinWorkDouble * region, int size, CoinWorkDouble value)
{
int i;
for (i = 0; i < size; i++)
region[i] = value;
}
void
multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
CoinWorkDouble * region2, CoinWorkDouble multiplier2)
{
int i;
if (multiplier1 == 1.0) {
if (multiplier2 == 1.0) {
for (i = 0; i < size; i++)
region2[i] = region1[i] + region2[i];
} else if (multiplier2 == -1.0) {
for (i = 0; i < size; i++)
region2[i] = region1[i] - region2[i];
} else if (multiplier2 == 0.0) {
for (i = 0; i < size; i++)
region2[i] = region1[i] ;
} else {
for (i = 0; i < size; i++)
region2[i] = region1[i] + multiplier2 * region2[i];
}
} else if (multiplier1 == -1.0) {
if (multiplier2 == 1.0) {
for (i = 0; i < size; i++)
region2[i] = -region1[i] + region2[i];
} else if (multiplier2 == -1.0) {
for (i = 0; i < size; i++)
region2[i] = -region1[i] - region2[i];
} else if (multiplier2 == 0.0) {
for (i = 0; i < size; i++)
region2[i] = -region1[i] ;
} else {
for (i = 0; i < size; i++)
region2[i] = -region1[i] + multiplier2 * region2[i];
}
} else if (multiplier1 == 0.0) {
if (multiplier2 == 1.0) {
// nothing to do
} else if (multiplier2 == -1.0) {
for (i = 0; i < size; i++)
region2[i] = -region2[i];
} else if (multiplier2 == 0.0) {
for (i = 0; i < size; i++)
region2[i] = 0.0;
} else {
for (i = 0; i < size; i++)
region2[i] = multiplier2 * region2[i];
}
} else {
if (multiplier2 == 1.0) {
for (i = 0; i < size; i++)
region2[i] = multiplier1 * region1[i] + region2[i];
} else if (multiplier2 == -1.0) {
for (i = 0; i < size; i++)
region2[i] = multiplier1 * region1[i] - region2[i];
} else if (multiplier2 == 0.0) {
for (i = 0; i < size; i++)
region2[i] = multiplier1 * region1[i] ;
} else {
for (i = 0; i < size; i++)
region2[i] = multiplier1 * region1[i] + multiplier2 * region2[i];
}
}
}
CoinWorkDouble
innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2)
{
int i;
CoinWorkDouble value = 0.0;
for (i = 0; i < size; i++)
value += region1[i] * region2[i];
return value;
}
void
getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2)
{
norm1 = 0.0;
norm2 = 0.0;
int i;
for (i = 0; i < size; i++) {
norm2 += region[i] * region[i];
norm1 = CoinMax(norm1, CoinAbs(region[i]));
}
}
#endif
#ifdef DEBUG_MEMORY
#include <malloc.h>
#include <stdio.h>
#include <stdlib.h>
typedef void (*NEW_HANDLER)();
static NEW_HANDLER new_handler; // function to call if `new' fails (cf. ARM p. 281)
// Allocate storage.
void *
operator new(size_t size)
{
void * p;
for (;;) {
p = malloc(size);
if (p) break; // success
else if (new_handler) new_handler(); // failure - try again (allow user to release some storage first)
else break; // failure - no retry
}
if (size > 1000000)
printf("Allocating memory of size %d\n", size);
return p;
}
// Deallocate storage.
void
operator delete(void *p)
{
free(p);
return;
}
void
operator delete [] (void *p)
{
free(p);
return;
}
#endif