source: trunk/Clp/src/ClpHelperFunctions.cpp @ 1367

Last change on this file since 1367 was 1367, checked in by forrest, 10 years ago

try and improve stability - also long double interior point

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.3 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/*
5    Note (JJF) I have added some operations on arrays even though they may
6    duplicate CoinDenseVector.  I think the use of templates was a mistake
7    as I don't think inline generic code can take as much advantage of
8    parallelism or machine architectures or memory hierarchies.
9
10*/
11#include <cfloat>
12#include <cstdlib>
13#include <cmath>
14#include "CoinHelperFunctions.hpp"
15#include "CoinFinite.hpp"
16double 
17maximumAbsElement(const double * region, int size)
18{
19  int i;
20  double maxValue=0.0;
21  for (i=0;i<size;i++) 
22    maxValue = CoinMax(maxValue,fabs(region[i]));
23  return maxValue;
24}
25void 
26setElements(double * region, int size, double value)
27{
28  int i;
29  for (i=0;i<size;i++) 
30    region[i]=value;
31}
32void 
33multiplyAdd(const double * region1, int size, double multiplier1,
34                 double * region2, double multiplier2)
35{
36  int i;
37  if (multiplier1==1.0) {
38    if (multiplier2==1.0) {
39      for (i=0;i<size;i++) 
40        region2[i] = region1[i] + region2[i];
41    } else if (multiplier2==-1.0) {
42      for (i=0;i<size;i++) 
43        region2[i] = region1[i] - region2[i];
44    } else if (multiplier2==0.0) {
45      for (i=0;i<size;i++) 
46        region2[i] = region1[i] ;
47    } else {
48      for (i=0;i<size;i++) 
49        region2[i] = region1[i] + multiplier2*region2[i];
50    }
51  } else if (multiplier1==-1.0) {
52    if (multiplier2==1.0) {
53      for (i=0;i<size;i++) 
54        region2[i] = -region1[i] + region2[i];
55    } else if (multiplier2==-1.0) {
56      for (i=0;i<size;i++) 
57        region2[i] = -region1[i] - region2[i];
58    } else if (multiplier2==0.0) {
59      for (i=0;i<size;i++) 
60        region2[i] = -region1[i] ;
61    } else {
62      for (i=0;i<size;i++) 
63        region2[i] = -region1[i] + multiplier2*region2[i];
64    }
65  } else if (multiplier1==0.0) {
66    if (multiplier2==1.0) {
67      // nothing to do
68    } else if (multiplier2==-1.0) {
69      for (i=0;i<size;i++) 
70        region2[i] =  -region2[i];
71    } else if (multiplier2==0.0) {
72      for (i=0;i<size;i++) 
73        region2[i] =  0.0;
74    } else {
75      for (i=0;i<size;i++) 
76        region2[i] =  multiplier2*region2[i];
77    }
78  } else {
79    if (multiplier2==1.0) {
80      for (i=0;i<size;i++) 
81        region2[i] = multiplier1*region1[i] + region2[i];
82    } else if (multiplier2==-1.0) {
83      for (i=0;i<size;i++) 
84        region2[i] = multiplier1*region1[i] - region2[i];
85    } else if (multiplier2==0.0) {
86      for (i=0;i<size;i++) 
87        region2[i] = multiplier1*region1[i] ;
88    } else {
89      for (i=0;i<size;i++) 
90        region2[i] = multiplier1*region1[i] + multiplier2*region2[i];
91    }
92  }
93}
94double 
95innerProduct(const double * region1, int size, const double * region2)
96{
97  int i;
98  double value=0.0;
99  for (i=0;i<size;i++)
100    value += region1[i]*region2[i];
101  return value;
102}
103void 
104getNorms(const double * region, int size, double & norm1, double & norm2)
105{
106  norm1 = 0.0;
107  norm2 = 0.0;
108  int i;
109  for (i=0;i<size;i++) {
110    norm2 += region[i]*region[i];
111    norm1 = CoinMax(norm1,fabs(region[i]));
112  }
113}
114#if COIN_LONG_WORK
115// For long double versions
116CoinWorkDouble
117maximumAbsElement(const CoinWorkDouble * region, int size)
118{
119  int i;
120  CoinWorkDouble maxValue=0.0;
121  for (i=0;i<size;i++) 
122    maxValue = CoinMax(maxValue,CoinAbs(region[i]));
123  return maxValue;
124}
125void 
126setElements(CoinWorkDouble * region, int size, CoinWorkDouble value)
127{
128  int i;
129  for (i=0;i<size;i++) 
130    region[i]=value;
131}
132void 
133multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
134                 CoinWorkDouble * region2, CoinWorkDouble multiplier2)
135{
136  int i;
137  if (multiplier1==1.0) {
138    if (multiplier2==1.0) {
139      for (i=0;i<size;i++) 
140        region2[i] = region1[i] + region2[i];
141    } else if (multiplier2==-1.0) {
142      for (i=0;i<size;i++) 
143        region2[i] = region1[i] - region2[i];
144    } else if (multiplier2==0.0) {
145      for (i=0;i<size;i++) 
146        region2[i] = region1[i] ;
147    } else {
148      for (i=0;i<size;i++) 
149        region2[i] = region1[i] + multiplier2*region2[i];
150    }
151  } else if (multiplier1==-1.0) {
152    if (multiplier2==1.0) {
153      for (i=0;i<size;i++) 
154        region2[i] = -region1[i] + region2[i];
155    } else if (multiplier2==-1.0) {
156      for (i=0;i<size;i++) 
157        region2[i] = -region1[i] - region2[i];
158    } else if (multiplier2==0.0) {
159      for (i=0;i<size;i++) 
160        region2[i] = -region1[i] ;
161    } else {
162      for (i=0;i<size;i++) 
163        region2[i] = -region1[i] + multiplier2*region2[i];
164    }
165  } else if (multiplier1==0.0) {
166    if (multiplier2==1.0) {
167      // nothing to do
168    } else if (multiplier2==-1.0) {
169      for (i=0;i<size;i++) 
170        region2[i] =  -region2[i];
171    } else if (multiplier2==0.0) {
172      for (i=0;i<size;i++) 
173        region2[i] =  0.0;
174    } else {
175      for (i=0;i<size;i++) 
176        region2[i] =  multiplier2*region2[i];
177    }
178  } else {
179    if (multiplier2==1.0) {
180      for (i=0;i<size;i++) 
181        region2[i] = multiplier1*region1[i] + region2[i];
182    } else if (multiplier2==-1.0) {
183      for (i=0;i<size;i++) 
184        region2[i] = multiplier1*region1[i] - region2[i];
185    } else if (multiplier2==0.0) {
186      for (i=0;i<size;i++) 
187        region2[i] = multiplier1*region1[i] ;
188    } else {
189      for (i=0;i<size;i++) 
190        region2[i] = multiplier1*region1[i] + multiplier2*region2[i];
191    }
192  }
193}
194CoinWorkDouble
195innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2)
196{
197  int i;
198  CoinWorkDouble value=0.0;
199  for (i=0;i<size;i++)
200    value += region1[i]*region2[i];
201  return value;
202}
203void 
204getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2)
205{
206  norm1 = 0.0;
207  norm2 = 0.0;
208  int i;
209  for (i=0;i<size;i++) {
210    norm2 += region[i]*region[i];
211    norm1 = CoinMax(norm1,CoinAbs(region[i]));
212  }
213}
214#endif
215#ifdef DEBUG_MEMORY
216#include <malloc.h>
217#include <stdio.h>
218#include <stdlib.h>
219
220typedef void (*NEW_HANDLER)();
221static NEW_HANDLER new_handler;                        // function to call if `new' fails (cf. ARM p. 281)
222
223// Allocate storage.
224void *
225operator new(size_t size)
226{
227  void * p;
228  for (;;)
229    {
230      p = malloc(size);
231      if      (p)           break;        // success
232      else if (new_handler) new_handler();   // failure - try again (allow user to release some storage first)
233      else                  break;        // failure - no retry
234    }
235  if (size>1000000)
236    printf("Allocating memory of size %d\n",size);
237  return p;
238}
239
240// Deallocate storage.
241void
242operator delete(void *p)
243{
244  free(p);
245  return;
246}
247void
248operator delete [] (void *p)
249{
250  free(p);
251  return;
252}
253#endif
Note: See TracBrowser for help on using the repository browser.