1/* $Id: ClpHelperFunctions.cpp 1753 2011-06-19 16:27:26Z stefan $ */
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 "CoinTypes.hpp"
18
19double
20maximumAbsElement(const double * region, int size)
21{
22 int i;
23 double maxValue = 0.0;
24 for (i = 0; i < size; i++)
25 maxValue = CoinMax(maxValue, fabs(region[i]));
26 return maxValue;
27}
28void
29setElements(double * region, int size, double value)
30{
31 int i;
32 for (i = 0; i < size; i++)
33 region[i] = value;
34}
35void
36multiplyAdd(const double * region1, int size, double multiplier1,
37 double * region2, double multiplier2)
38{
39 int i;
40 if (multiplier1 == 1.0) {
41 if (multiplier2 == 1.0) {
42 for (i = 0; i < size; i++)
43 region2[i] = region1[i] + region2[i];
44 } else if (multiplier2 == -1.0) {
45 for (i = 0; i < size; i++)
46 region2[i] = region1[i] - region2[i];
47 } else if (multiplier2 == 0.0) {
48 for (i = 0; i < size; i++)
49 region2[i] = region1[i] ;
50 } else {
51 for (i = 0; i < size; i++)
52 region2[i] = region1[i] + multiplier2 * region2[i];
53 }
54 } else if (multiplier1 == -1.0) {
55 if (multiplier2 == 1.0) {
56 for (i = 0; i < size; i++)
57 region2[i] = -region1[i] + region2[i];
58 } else if (multiplier2 == -1.0) {
59 for (i = 0; i < size; i++)
60 region2[i] = -region1[i] - region2[i];
61 } else if (multiplier2 == 0.0) {
62 for (i = 0; i < size; i++)
63 region2[i] = -region1[i] ;
64 } else {
65 for (i = 0; i < size; i++)
66 region2[i] = -region1[i] + multiplier2 * region2[i];
67 }
68 } else if (multiplier1 == 0.0) {
69 if (multiplier2 == 1.0) {
70 // nothing to do
71 } else if (multiplier2 == -1.0) {
72 for (i = 0; i < size; i++)
73 region2[i] = -region2[i];
74 } else if (multiplier2 == 0.0) {
75 for (i = 0; i < size; i++)
76 region2[i] = 0.0;
77 } else {
78 for (i = 0; i < size; i++)
79 region2[i] = multiplier2 * region2[i];
80 }
81 } else {
82 if (multiplier2 == 1.0) {
83 for (i = 0; i < size; i++)
84 region2[i] = multiplier1 * region1[i] + region2[i];
85 } else if (multiplier2 == -1.0) {
86 for (i = 0; i < size; i++)
87 region2[i] = multiplier1 * region1[i] - region2[i];
88 } else if (multiplier2 == 0.0) {
89 for (i = 0; i < size; i++)
90 region2[i] = multiplier1 * region1[i] ;
91 } else {
92 for (i = 0; i < size; i++)
93 region2[i] = multiplier1 * region1[i] + multiplier2 * region2[i];
94 }
95 }
96}
97double
98innerProduct(const double * region1, int size, const double * region2)
99{
100 int i;
101 double value = 0.0;
102 for (i = 0; i < size; i++)
103 value += region1[i] * region2[i];
104 return value;
105}
106void
107getNorms(const double * region, int size, double & norm1, double & norm2)
108{
109 norm1 = 0.0;
110 norm2 = 0.0;
111 int i;
112 for (i = 0; i < size; i++) {
113 norm2 += region[i] * region[i];
114 norm1 = CoinMax(norm1, fabs(region[i]));
115 }
116}
117#if COIN_LONG_WORK
118// For long double versions
119CoinWorkDouble
120maximumAbsElement(const CoinWorkDouble * region, int size)
121{
122 int i;
123 CoinWorkDouble maxValue = 0.0;
124 for (i = 0; i < size; i++)
125 maxValue = CoinMax(maxValue, CoinAbs(region[i]));
126 return maxValue;
127}
128void
129setElements(CoinWorkDouble * region, int size, CoinWorkDouble value)
130{
131 int i;
132 for (i = 0; i < size; i++)
133 region[i] = value;
134}
135void
136multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
137 CoinWorkDouble * region2, CoinWorkDouble multiplier2)
138{
139 int i;
140 if (multiplier1 == 1.0) {
141 if (multiplier2 == 1.0) {
142 for (i = 0; i < size; i++)
143 region2[i] = region1[i] + region2[i];
144 } else if (multiplier2 == -1.0) {
145 for (i = 0; i < size; i++)
146 region2[i] = region1[i] - region2[i];
147 } else if (multiplier2 == 0.0) {
148 for (i = 0; i < size; i++)
149 region2[i] = region1[i] ;
150 } else {
151 for (i = 0; i < size; i++)
152 region2[i] = region1[i] + multiplier2 * region2[i];
153 }
154 } else if (multiplier1 == -1.0) {
155 if (multiplier2 == 1.0) {
156 for (i = 0; i < size; i++)
157 region2[i] = -region1[i] + region2[i];
158 } else if (multiplier2 == -1.0) {
159 for (i = 0; i < size; i++)
160 region2[i] = -region1[i] - region2[i];
161 } else if (multiplier2 == 0.0) {
162 for (i = 0; i < size; i++)
163 region2[i] = -region1[i] ;
164 } else {
165 for (i = 0; i < size; i++)
166 region2[i] = -region1[i] + multiplier2 * region2[i];
167 }
168 } else if (multiplier1 == 0.0) {
169 if (multiplier2 == 1.0) {
170 // nothing to do
171 } else if (multiplier2 == -1.0) {
172 for (i = 0; i < size; i++)
173 region2[i] = -region2[i];
174 } else if (multiplier2 == 0.0) {
175 for (i = 0; i < size; i++)
176 region2[i] = 0.0;
177 } else {
178 for (i = 0; i < size; i++)
179 region2[i] = multiplier2 * region2[i];
180 }
181 } else {
182 if (multiplier2 == 1.0) {
183 for (i = 0; i < size; i++)
184 region2[i] = multiplier1 * region1[i] + region2[i];
185 } else if (multiplier2 == -1.0) {
186 for (i = 0; i < size; i++)
187 region2[i] = multiplier1 * region1[i] - region2[i];
188 } else if (multiplier2 == 0.0) {
189 for (i = 0; i < size; i++)
190 region2[i] = multiplier1 * region1[i] ;
191 } else {
192 for (i = 0; i < size; i++)
193 region2[i] = multiplier1 * region1[i] + multiplier2 * region2[i];
194 }
195 }
196}
197CoinWorkDouble
198innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2)
199{
200 int i;
201 CoinWorkDouble value = 0.0;
202 for (i = 0; i < size; i++)
203 value += region1[i] * region2[i];
204 return value;
205}
206void
207getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2)
208{
209 norm1 = 0.0;
210 norm2 = 0.0;
211 int i;
212 for (i = 0; i < size; i++) {
213 norm2 += region[i] * region[i];
214 norm1 = CoinMax(norm1, CoinAbs(region[i]));
215 }
216}
217#endif
218#ifdef DEBUG_MEMORY
219#include <malloc.h>
220#include <stdio.h>
221#include <stdlib.h>
222
223typedef void (*NEW_HANDLER)();
224static NEW_HANDLER new_handler; // function to call if `new' fails (cf. ARM p. 281)
225
226// Allocate storage.
227void *
228operator new(size_t size)
229{
230 void * p;
231 for (;;) {
232 p = malloc(size);
233 if (p) break; // success
234 else if (new_handler) new_handler(); // failure - try again (allow user to release some storage first)
235 else break; // failure - no retry
236 }
237 if (size > 1000000)
238 printf("Allocating memory of size %d\n", size);
239 return p;
240}
241
242// Deallocate storage.
243void
244operator delete(void *p)
245{
246 free(p);
247 return;
248}
249void
250operator delete [] (void *p)
251{
252 free(p);
253 return;
254}
255#endif
256