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

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

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