source: trunk/Clp/src/ClpMain.cpp @ 2449

Last change on this file since 2449 was 2426, checked in by stefan, 9 months ago

replace remaining cilk_grainsize by cilk grainsize, closes #68

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 10.8 KB
Line 
1/* $Id: ClpMain.cpp 2426 2019-03-11 15:13:16Z stefan $ */
2// Copyright (C) 2002, 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#include "CoinPragma.hpp"
7
8#include "AbcCommon.hpp"
9#include "ClpSimplex.hpp"
10#ifdef ABC_INHERIT
11#include "AbcSimplex.hpp"
12#endif
13#ifndef ABC_INHERIT
14void ClpMain0(ClpSimplex *models);
15int ClpMain1(int argc, const char *argv[], ClpSimplex *model);
16#else
17void ClpMain0(AbcSimplex *models);
18int ClpMain1(int argc, const char *argv[], AbcSimplex *model);
19#endif
20//#define CILK_TEST
21#ifdef CILK_TEST
22static void cilkTest();
23#endif
24//#define LAPACK_TEST
25/*
26  Somehow with some BLAS we get multithreaded by default
27  For 99.99% of problems this is not a good idea.
28  The openblas_set_num_threads(1) seems to work even with other blas
29 */
30#if CLP_USE_OPENBLAS
31extern "C" {
32void openblas_set_num_threads(int num_threads);
33}
34#endif
35#ifdef LAPACK_TEST
36#include "/include/lapacke.h"
37#ifndef COIN_FACTORIZATION_DENSE_CODE
38#define COIN_FACTORIZATION_DENSE_CODE 1
39#endif
40#if COIN_FACTORIZATION_DENSE_CODE == 1
41// using simple lapack interface
42extern "C" {
43void openblas_set_num_threads(int num_threads);
44#if 0
45  /** LAPACK Fortran subroutine DGETRF. */
46  void LAPACK_dgetrf(int * m, int *n,
47                               double *A, int *ldA,
48                               int * ipiv, int *info);
49  /** LAPACK Fortran subroutine DGETRS. */
50  void LAPACK_dgetrs(char *trans, int *n,
51                               int *nrhs, const double *A, int *ldA,
52                     int * ipiv, double *B, int *ldB, int *info);
53  //    LAPACK_dgetrf(&N, &N, m, &LDA,ipiv, &info);
54#endif
55}
56int test_lapack(int n)
57{
58  int *ipiv;
59  int info;
60  int i, j;
61  double *m, *x, *y;
62
63  int LDB, LDA, N, NRHS;
64  char transp = 'N';
65
66  m = (double *)malloc(sizeof(double) * n * n);
67  x = (double *)malloc(sizeof(double) * n);
68  y = (double *)malloc(sizeof(double) * n);
69  ipiv = (int *)malloc(sizeof(int) * n);
70
71  for (i = 0; i < n; ++i) {
72    x[i] = 1.0;
73    for (j = 0; j < n; ++j) {
74      m[i * n + j] = (rand() % 100 + 1) / 10.0;
75      //      printf("m[%d,%d]=%lf\n",i,j, m[i*n+j]);
76    }
77  }
78
79  /* test cblas.h */
80  //cblas_dgemv(CblasColMajor, CblasNoTrans, n, n, 1.0, m, n,
81  //          x, 1, 0.0, y, 1);
82
83  //  for (i=0; i<n; ++i)  printf("x[%d]=%lf\n",i, x[i]);
84  //for (i=0; i<n; ++i)  printf("y[%d]=%lf\n",i, y[i]);
85
86  LDB = n;
87  LDA = n;
88  N = n;
89  NRHS = 1;
90  info = 0;
91
92  LAPACK_dgetrf(&N, &N, m, &LDA, ipiv, &info);
93
94  if (info != 0)
95    fprintf(stderr, "dgetrf failure with error %d\n", info);
96
97  LAPACK_dgetrs(&transp, &N, &NRHS, m, &LDA, ipiv, y, &LDB, &info);
98
99  if (info != 0)
100    fprintf(stderr, "failure with error %d\n", info);
101  //  for (i=0; i<n; ++i) printf("%lf\n", y[i]);
102
103  free(m);
104  free(x);
105  free(y);
106  free(ipiv);
107  return 0;
108}
109#elif COIN_FACTORIZATION_DENSE_CODE == 2
110// C interface
111enum CBLAS_ORDER { CblasRowMajor = 101,
112  CblasColMajor = 102 };
113enum CBLAS_TRANSPOSE { CblasNoTrans = 111,
114  CblasTrans = 112 };
115extern "C" {
116int clapack_dgetrf(const enum CBLAS_ORDER Order, const int M, const int N, double *A, const int lda, int *ipiv);
117int clapack_dgetrs(const enum CBLAS_ORDER Order,
118  const enum CBLAS_TRANSPOSE Trans,
119  const int N, const int NRHS,
120  const double *A, const int lda, const int *ipiv, double *B,
121  const int ldb);
122}
123int test_lapack(int n)
124{
125  int *ipiv;
126  int info;
127  int i, j;
128  double *m, *x, *y;
129
130  int LDB, LDA, N, NRHS;
131  char transp = 'N';
132
133  m = (double *)malloc(sizeof(double) * n * n);
134  x = (double *)malloc(sizeof(double) * n);
135  y = (double *)malloc(sizeof(double) * n);
136  ipiv = (int *)malloc(sizeof(int) * n);
137
138  for (i = 0; i < n; ++i) {
139    x[i] = 1.0;
140    for (j = 0; j < n; ++j) {
141      m[i * n + j] = (rand() % 100 + 1) / 10.0;
142      //      printf("m[%d,%d]=%lf\n",i,j, m[i*n+j]);
143    }
144  }
145
146  /* test cblas.h */
147  //cblas_dgemv(CblasColMajor, CblasNoTrans, n, n, 1.0, m, n,
148  //          x, 1, 0.0, y, 1);
149
150  //  for (i=0; i<n; ++i)  printf("x[%d]=%lf\n",i, x[i]);
151  //for (i=0; i<n; ++i)  printf("y[%d]=%lf\n",i, y[i]);
152
153  LDB = n;
154  LDA = n;
155  N = n;
156  NRHS = 1;
157  info = clapack_dgetrf(CblasColMajor, n, n, m, n, ipiv);
158
159  if (info != 0)
160    fprintf(stderr, "dgetrf failure with error %d\n", info);
161
162  clapack_dgetrs(CblasColMajor, CblasNoTrans, n, 1, m, n, ipiv, y, n);
163
164  if (info != 0)
165    fprintf(stderr, "failure with error %d\n", info);
166  //  for (i=0; i<n; ++i) printf("%lf\n", y[i]);
167
168  free(m);
169  free(x);
170  free(y);
171  free(ipiv);
172  return 0;
173}
174#endif
175#endif
176//#define CLP_MALLOC_STATISTICS
177#ifdef CLP_MALLOC_STATISTICS
178#include <malloc.h>
179#include <exception>
180#include <new>
181static double malloc_times = 0.0;
182static double malloc_total = 0.0;
183static int malloc_amount[] = { 0, 32, 128, 256, 1024, 4096, 16384, 65536, 262144, COIN_INT_MAX };
184static int malloc_n = 10;
185double malloc_counts[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
186bool malloc_counts_on = true;
187void *operator new(size_t size) throw(std::bad_alloc)
188{
189  malloc_times++;
190  malloc_total += size;
191  int i;
192  for (i = 0; i < malloc_n; i++) {
193    if ((int)size <= malloc_amount[i]) {
194      malloc_counts[i]++;
195      break;
196    }
197  }
198  if (size > 100000) {
199    printf("allocating %ld bytes\n", size);
200  }
201  void *p = malloc(size);
202  return p;
203}
204void operator delete(void *p) throw()
205{
206  free(p);
207}
208static void malloc_stats2()
209{
210  double average = malloc_total / malloc_times;
211  printf("count %g bytes %g - average %g\n", malloc_times, malloc_total, average);
212  for (int i = 0; i < malloc_n; i++)
213    printf("%g ", malloc_counts[i]);
214  printf("\n");
215  malloc_times = 0.0;
216  malloc_total = 0.0;
217  memset(malloc_counts, 0, sizeof(malloc_counts));
218  // print results
219}
220#endif //CLP_MALLOC_STATISTICS
221int
222#if defined(_MSC_VER)
223  __cdecl
224#endif // _MSC_VER
225  main(int argc, const char *argv[])
226{
227#ifdef CILK_TEST
228  cilkTest();
229#endif
230#if CLP_USE_OPENBLAS
231  openblas_set_num_threads(CLP_USE_OPENBLAS);
232#endif
233#ifdef LAPACK_TEST
234  //void openblas_set_num_threads(int num_threads);
235  openblas_set_num_threads(1);
236  if (argc < 2) {
237    printf("Error - need size of matrix for lapack test\n");
238    return 1;
239  }
240  int n = atoi(argv[1]);
241  printf("n=%d\n", n);
242  if (argc > 2) {
243    int nThreads = atoi(argv[2]);
244    printf("Using %d threads\n", nThreads);
245    openblas_set_num_threads(nThreads);
246  }
247  test_lapack(n);
248  return 0;
249#endif
250#ifndef ABC_INHERIT
251  ClpSimplex *models = new ClpSimplex[1];
252#else
253  AbcSimplex *models = new AbcSimplex[1];
254#endif
255  std::cout << "Coin LP version " << CLP_VERSION
256            << ", build " << __DATE__ << std::endl;
257  // Print command line
258  if (argc > 1) {
259    printf("command line - ");
260    for (int i = 0; i < argc; i++)
261      printf("%s ", argv[i]);
262    printf("\n");
263  }
264  ClpMain0(models);
265  int returnCode = ClpMain1(argc, argv, models);
266  delete[] models;
267  return returnCode;
268}
269/*
270  Version 1.00.00 October 13 2004.
271  1.00.01 October 18.  Added basis handling helped/prodded by Thorsten Koch.
272  Also modifications to make faster with sbb (I hope I haven't broken anything).
273  1.00.02 March 21 2005.  Redid ClpNonLinearCost to save memory also redid
274  createRim to try and improve cache characteristics.
275  1.00.03 April 8 2005.  Added Volume algorithm as crash and made code more
276  robust on testing.  Also added "either" and "tune" algorithm.
277  1.01.01 April 12 2005.  Decided to go to different numbering.  Backups will
278  be last 2 digits while middle 2 are for improvements.  Still take a long
279  time to get to 2.00.01
280  1.01.02 May 4 2005.  Will be putting in many changes - so saving stable version
281  1.02.01 May 6 2005.  Lots of changes to try and make faster and more stable in
282  branch and cut.
283  1.02.02 May 19 2005.  Stuff for strong branching and some improvements to simplex
284  1.03.01 May 24 2006.  Lots done but I can't remember what!
285  1.03.03 June 13 2006.  For clean up after dual perturbation
286  1.04.01 June 26 2007.  Lots of changes but I got lazy
287  1.05.00 June 27 2007.  This is trunk so when gets to stable will be 1.5
288  1.11.00 November 5 2009 (Guy Fawkes) - OSL factorization and better ordering
289 */
290#ifdef CILK_TEST
291// -*- C++ -*-
292
293/*
294 * cilk-for.cilk
295 *
296 * Copyright (c) 2007-2008 Cilk Arts, Inc.  55 Cambridge Street,
297 * Burlington, MA 01803.  Patents pending.  All rights reserved. You may
298 * freely use the sample code to guide development of your own works,
299 * provided that you reproduce this notice in any works you make that
300 * use the sample code.  This sample code is provided "AS IS" without
301 * warranty of any kind, either express or implied, including but not
302 * limited to any implied warranty of non-infringement, merchantability
303 * or fitness for a particular purpose.  In no event shall Cilk Arts,
304 * Inc. be liable for any direct, indirect, special, or consequential
305 * damages, or any other damages whatsoever, for any use of or reliance
306 * on this sample code, including, without limitation, any lost
307 * opportunity, lost profits, business interruption, loss of programs or
308 * data, even if expressly advised of or otherwise aware of the
309 * possibility of such damages, whether in an action of contract,
310 * negligence, tort, or otherwise.
311 *
312 * This file demonstrates a Cilk++ for loop
313 */
314
315#include <cilk/cilk.h>
316//#include <cilk/cilkview.h>
317#include <cilk/reducer_max.h>
318#include <cstdlib>
319#include <iostream>
320
321// cilk_for granularity.
322#define CILK_FOR_GRAINSIZE 128
323
324double dowork(double i)
325{
326  // Waste time:
327  int j;
328  double k = i;
329  for (j = 0; j < 50000; ++j) {
330    k += k / ((j + 1) * (k + 1));
331  }
332
333  return k;
334}
335static void doSomeWork(double *a, int low, int high)
336{
337  if (high - low > 300) {
338    int mid = (high + low) >> 1;
339    cilk_spawn doSomeWork(a, low, mid);
340    doSomeWork(a, mid, high);
341    cilk_sync;
342  } else {
343    for (int i = low; i < high; ++i) {
344      a[i] = dowork(a[i]);
345    }
346  }
347}
348
349void cilkTest()
350{
351  unsigned int n = 10000;
352  //cilk::cilkview cv;
353
354  double *a = new double[n];
355
356  for (unsigned int i = 0; i < n; i++) {
357    // Populate A
358    a[i] = (double)((i * i) % 1024 + 512) / 512;
359  }
360
361  std::cout << "Iterating over " << n << " integers" << std::endl;
362
363  //cv.start();
364#if 1
365  //#pragma cilk grainsize = CILK_FOR_GRAINSIZE
366  cilk_for(unsigned int i = 0; i < n; ++i)
367  {
368    a[i] = dowork(a[i]);
369  }
370#else
371  doSomeWork(a, 0, n);
372#endif
373  int *which = new int[n];
374  unsigned int n2 = n >> 1;
375  for (int i = 0; i < n2; i++)
376    which[i] = n - 2 * i;
377  cilk::reducer_max_index< int, double > maximumIndex(-1, 0.0);
378  cilk_for(unsigned int i = 0; i < n2; ++i)
379  {
380    int iWhich = which[i];
381    maximumIndex.calc_max(iWhich, a[iWhich]);
382  }
383  int bestIndex = maximumIndex.get_index();
384  int bestIndex2 = -1;
385  double largest = 0.0;
386  cilk_for(unsigned int i = 0; i < n2; ++i)
387  {
388    int iWhich = which[i];
389    if (a[iWhich] > largest) {
390      bestIndex2 = iWhich;
391      largest = a[iWhich];
392    }
393  }
394  assert(bestIndex == bestIndex2);
395  //cv.stop();
396  //cv.dump("cilk-for-results", false);
397
398  //std::cout << cv.accumulated_milliseconds() / 1000.f << " seconds" << std::endl;
399
400  exit(0);
401}
402#endif
403
404/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
405*/
Note: See TracBrowser for help on using the repository browser.