/farmR/src/glpbfd.c

https://code.google.com/p/javawfm/ · C · 419 lines · 200 code · 17 blank · 202 comment · 66 complexity · b6b22c0568fc7067cc5511eeda83e93b MD5 · raw file

  1. /* glpbfd.c (LP basis factorization driver) */
  2. /***********************************************************************
  3. * This code is part of GLPK (GNU Linear Programming Kit).
  4. *
  5. * Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
  6. * Department for Applied Informatics, Moscow Aviation Institute,
  7. * Moscow, Russia. All rights reserved. E-mail: <mao@mai2.rcnet.ru>.
  8. *
  9. * GLPK is free software: you can redistribute it and/or modify it
  10. * under the terms of the GNU General Public License as published by
  11. * the Free Software Foundation, either version 3 of the License, or
  12. * (at your option) any later version.
  13. *
  14. * GLPK is distributed in the hope that it will be useful, but WITHOUT
  15. * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16. * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
  17. * License for more details.
  18. *
  19. * You should have received a copy of the GNU General Public License
  20. * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
  21. ***********************************************************************/
  22. #define _GLPBFD_PRIVATE
  23. #include "glpbfd.h"
  24. #include "glplib.h"
  25. #define xfault xerror
  26. /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
  27. #define M_MAX 100000000 /* = 100*10^6 */
  28. /* maximal order of the basis matrix */
  29. /***********************************************************************
  30. * NAME
  31. *
  32. * bfd_create_it - create LP basis factorization
  33. *
  34. * SYNOPSIS
  35. *
  36. * #include "glpbfd.h"
  37. * BFD *bfd_create_it(void);
  38. *
  39. * DESCRIPTION
  40. *
  41. * The routine bfd_create_it creates a program object, which represents
  42. * a factorization of LP basis.
  43. *
  44. * RETURNS
  45. *
  46. * The routine bfd_create_it returns a pointer to the object created. */
  47. BFD *bfd_create_it(void)
  48. { BFD *bfd;
  49. bfd = xmalloc(sizeof(BFD));
  50. bfd->valid = 0;
  51. bfd->type = BFD_TFT;
  52. bfd->fhv = NULL;
  53. bfd->lpf = NULL;
  54. bfd->lu_size = 0;
  55. bfd->piv_tol = 0.10;
  56. bfd->piv_lim = 4;
  57. bfd->suhl = 1;
  58. bfd->eps_tol = 1e-15;
  59. bfd->max_gro = 1e+10;
  60. bfd->nfs_max = 50;
  61. bfd->upd_tol = 1e-6;
  62. bfd->nrs_max = 50;
  63. bfd->rs_size = 1000;
  64. bfd->upd_lim = -1;
  65. bfd->upd_cnt = 0;
  66. return bfd;
  67. }
  68. /***********************************************************************
  69. * NAME
  70. *
  71. * bfd_factorize - compute LP basis factorization
  72. *
  73. * SYNOPSIS
  74. *
  75. * #include "glpbfd.h"
  76. * int bfd_factorize(BFD *bfd, int m, int bh[], int (*col)(void *info,
  77. * int j, int ind[], double val[]), void *info);
  78. *
  79. * DESCRIPTION
  80. *
  81. * The routine bfd_factorize computes the factorization of the basis
  82. * matrix B specified by the routine col.
  83. *
  84. * The parameter bfd specified the basis factorization data structure
  85. * created with the routine bfd_create_it.
  86. *
  87. * The parameter m specifies the order of B, m > 0.
  88. *
  89. * The array bh specifies the basis header: bh[j], 1 <= j <= m, is the
  90. * number of j-th column of B in some original matrix. The array bh is
  91. * optional and can be specified as NULL.
  92. *
  93. * The formal routine col specifies the matrix B to be factorized. To
  94. * obtain j-th column of A the routine bfd_factorize calls the routine
  95. * col with the parameter j (1 <= j <= n). In response the routine col
  96. * should store row indices and numerical values of non-zero elements
  97. * of j-th column of B to locations ind[1,...,len] and val[1,...,len],
  98. * respectively, where len is the number of non-zeros in j-th column
  99. * returned on exit. Neither zero nor duplicate elements are allowed.
  100. *
  101. * The parameter info is a transit pointer passed to the routine col.
  102. *
  103. * RETURNS
  104. *
  105. * 0 The factorization has been successfully computed.
  106. *
  107. * BFD_ESING
  108. * The specified matrix is singular within the working precision.
  109. *
  110. * BFD_ECOND
  111. * The specified matrix is ill-conditioned.
  112. *
  113. * For more details see comments to the routine luf_factorize. */
  114. int bfd_factorize(BFD *bfd, int m, const int bh[], int (*col)
  115. (void *info, int j, int ind[], double val[]), void *info)
  116. { LUF *luf;
  117. int nov, ret;
  118. if (m < 1)
  119. xfault("bfd_factorize: m = %d; invalid parameter\n", m);
  120. if (m > M_MAX)
  121. xfault("bfd_factorize: m = %d; matrix too big\n", m);
  122. /* invalidate the factorization */
  123. bfd->valid = 0;
  124. /* create the factorization, if necessary */
  125. nov = 0;
  126. switch (bfd->type)
  127. { case BFD_TFT:
  128. if (bfd->lpf != NULL)
  129. lpf_delete_it(bfd->lpf), bfd->lpf = NULL;
  130. if (bfd->fhv == NULL)
  131. bfd->fhv = fhv_create_it(), nov = 1;
  132. break;
  133. case BFD_TBG:
  134. case BFD_TGR:
  135. if (bfd->fhv != NULL)
  136. fhv_delete_it(bfd->fhv), bfd->fhv = NULL;
  137. if (bfd->lpf == NULL)
  138. bfd->lpf = lpf_create_it(), nov = 1;
  139. break;
  140. default:
  141. xfault("bfd_factorize: bfd->type = %d; invalid factorizatio"
  142. "n type\n", bfd->type);
  143. }
  144. /* set control parameters specific to LUF */
  145. if (bfd->fhv != NULL)
  146. luf = bfd->fhv->luf;
  147. else if (bfd->lpf != NULL)
  148. luf = bfd->lpf->luf;
  149. else
  150. xassert(bfd != bfd);
  151. if (nov) luf->new_sva = bfd->lu_size;
  152. luf->piv_tol = bfd->piv_tol;
  153. luf->piv_lim = bfd->piv_lim;
  154. luf->suhl = bfd->suhl;
  155. luf->eps_tol = bfd->eps_tol;
  156. luf->max_gro = bfd->max_gro;
  157. /* set control parameters specific to FHV */
  158. if (bfd->fhv != NULL)
  159. { if (nov) bfd->fhv->hh_max = bfd->nfs_max;
  160. bfd->fhv->upd_tol = bfd->upd_tol;
  161. }
  162. /* set control parameters specific to LPF */
  163. if (bfd->lpf != NULL)
  164. { if (nov) bfd->lpf->n_max = bfd->nrs_max;
  165. if (nov) bfd->lpf->v_size = bfd->rs_size;
  166. }
  167. /* try to factorize the basis matrix */
  168. if (bfd->fhv != NULL)
  169. { switch (fhv_factorize(bfd->fhv, m, col, info))
  170. { case 0:
  171. break;
  172. case FHV_ESING:
  173. ret = BFD_ESING;
  174. goto done;
  175. case FHV_ECOND:
  176. ret = BFD_ECOND;
  177. goto done;
  178. default:
  179. xassert(bfd != bfd);
  180. }
  181. }
  182. else if (bfd->lpf != NULL)
  183. { switch (lpf_factorize(bfd->lpf, m, bh, col, info))
  184. { case 0:
  185. /* set the Schur complement update type */
  186. switch (bfd->type)
  187. { case BFD_TBG:
  188. /* Bartels-Golub update */
  189. bfd->lpf->scf->t_opt = SCF_TBG;
  190. break;
  191. case BFD_TGR:
  192. /* Givens rotation update */
  193. bfd->lpf->scf->t_opt = SCF_TGR;
  194. break;
  195. default:
  196. xassert(bfd != bfd);
  197. }
  198. break;
  199. case LPF_ESING:
  200. ret = BFD_ESING;
  201. goto done;
  202. case LPF_ECOND:
  203. ret = BFD_ECOND;
  204. goto done;
  205. default:
  206. xassert(bfd != bfd);
  207. }
  208. }
  209. else
  210. xassert(bfd != bfd);
  211. /* the basis matrix has been successfully factorized */
  212. bfd->valid = 1;
  213. bfd->upd_cnt = 0;
  214. ret = 0;
  215. done: /* return to the calling program */
  216. return ret;
  217. }
  218. /***********************************************************************
  219. * NAME
  220. *
  221. * bfd_ftran - perform forward transformation (solve system B*x = b)
  222. *
  223. * SYNOPSIS
  224. *
  225. * #include "glpbfd.h"
  226. * void bfd_ftran(BFD *bfd, double x[]);
  227. *
  228. * DESCRIPTION
  229. *
  230. * The routine bfd_ftran performs forward transformation, i.e. solves
  231. * the system B*x = b, where B is the basis matrix, x is the vector of
  232. * unknowns to be computed, b is the vector of right-hand sides.
  233. *
  234. * On entry elements of the vector b should be stored in dense format
  235. * in locations x[1], ..., x[m], where m is the number of rows. On exit
  236. * the routine stores elements of the vector x in the same locations. */
  237. void bfd_ftran(BFD *bfd, double x[])
  238. { if (!bfd->valid)
  239. xfault("bfd_ftran: the factorization is not valid\n");
  240. if (bfd->fhv != NULL)
  241. fhv_ftran(bfd->fhv, x);
  242. else if (bfd->lpf != NULL)
  243. lpf_ftran(bfd->lpf, x);
  244. else
  245. xassert(bfd != bfd);
  246. return;
  247. }
  248. /***********************************************************************
  249. * NAME
  250. *
  251. * bfd_btran - perform backward transformation (solve system B'*x = b)
  252. *
  253. * SYNOPSIS
  254. *
  255. * #include "glpbfd.h"
  256. * void bfd_btran(BFD *bfd, double x[]);
  257. *
  258. * DESCRIPTION
  259. *
  260. * The routine bfd_btran performs backward transformation, i.e. solves
  261. * the system B'*x = b, where B' is a matrix transposed to the basis
  262. * matrix B, x is the vector of unknowns to be computed, b is the vector
  263. * of right-hand sides.
  264. *
  265. * On entry elements of the vector b should be stored in dense format
  266. * in locations x[1], ..., x[m], where m is the number of rows. On exit
  267. * the routine stores elements of the vector x in the same locations. */
  268. void bfd_btran(BFD *bfd, double x[])
  269. { if (!bfd->valid)
  270. xfault("bfd_btran: the factorization is not valid\n");
  271. if (bfd->fhv != NULL)
  272. fhv_btran(bfd->fhv, x);
  273. else if (bfd->lpf != NULL)
  274. lpf_btran(bfd->lpf, x);
  275. else
  276. xassert(bfd != bfd);
  277. return;
  278. }
  279. /***********************************************************************
  280. * NAME
  281. *
  282. * bfd_update_it - update LP basis factorization
  283. *
  284. * SYNOPSIS
  285. *
  286. * #include "glpbfd.h"
  287. * int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[],
  288. * const double val[]);
  289. *
  290. * DESCRIPTION
  291. *
  292. * The routine bfd_update_it updates the factorization of the basis
  293. * matrix B after replacing its j-th column by a new vector.
  294. *
  295. * The parameter j specifies the number of column of B, which has been
  296. * replaced, 1 <= j <= m, where m is the order of B.
  297. *
  298. * The parameter bh specifies the basis header entry for the new column
  299. * of B, which is the number of the new column in some original matrix.
  300. * This parameter is optional and can be specified as 0.
  301. *
  302. * Row indices and numerical values of non-zero elements of the new
  303. * column of B should be placed in locations ind[1], ..., ind[len] and
  304. * val[1], ..., val[len], resp., where len is the number of non-zeros
  305. * in the column. Neither zero nor duplicate elements are allowed.
  306. *
  307. * RETURNS
  308. *
  309. * 0 The factorization has been successfully updated.
  310. *
  311. * BFD_ESING
  312. * New basis matrix is singular within the working precision.
  313. *
  314. * BFD_ECHECK
  315. * The factorization is inaccurate.
  316. *
  317. * BFD_ELIMIT
  318. * Factorization update limit has been reached.
  319. *
  320. * BFD_EROOM
  321. * Overflow of the sparse vector area.
  322. *
  323. * In case of non-zero return code the factorization becomes invalid.
  324. * It should not be used until it has been recomputed with the routine
  325. * bfd_factorize. */
  326. int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[],
  327. const double val[])
  328. { int ret;
  329. if (!bfd->valid)
  330. xfault("bfd_update_it: the factorization is not valid\n");
  331. /* try to update the factorization */
  332. if (bfd->fhv != NULL)
  333. { switch (fhv_update_it(bfd->fhv, j, len, ind, val))
  334. { case 0:
  335. break;
  336. case FHV_ESING:
  337. bfd->valid = 0;
  338. ret = BFD_ESING;
  339. goto done;
  340. case FHV_ECHECK:
  341. bfd->valid = 0;
  342. ret = BFD_ECHECK;
  343. goto done;
  344. case FHV_ELIMIT:
  345. bfd->valid = 0;
  346. ret = BFD_ELIMIT;
  347. goto done;
  348. case FHV_EROOM:
  349. bfd->valid = 0;
  350. ret = BFD_EROOM;
  351. goto done;
  352. default:
  353. xassert(bfd != bfd);
  354. }
  355. }
  356. else if (bfd->lpf != NULL)
  357. { switch (lpf_update_it(bfd->lpf, j, bh, len, ind, val))
  358. { case 0:
  359. break;
  360. case LPF_ESING:
  361. bfd->valid = 0;
  362. ret = BFD_ESING;
  363. goto done;
  364. case LPF_ELIMIT:
  365. bfd->valid = 0;
  366. ret = BFD_ELIMIT;
  367. goto done;
  368. default:
  369. xassert(bfd != bfd);
  370. }
  371. }
  372. else
  373. xassert(bfd != bfd);
  374. /* the factorization has been successfully updated */
  375. /* increase the update count */
  376. bfd->upd_cnt++;
  377. ret = 0;
  378. done: /* return to the calling program */
  379. return ret;
  380. }
  381. /***********************************************************************
  382. * NAME
  383. *
  384. * bfd_delete_it - delete LP basis factorization
  385. *
  386. * SYNOPSIS
  387. *
  388. * #include "glpbfd.h"
  389. * void bfd_delete_it(BFD *bfd);
  390. *
  391. * DESCRIPTION
  392. *
  393. * The routine bfd_delete_it deletes LP basis factorization specified
  394. * by the parameter fhv and frees all memory allocated to this program
  395. * object. */
  396. void bfd_delete_it(BFD *bfd)
  397. { if (bfd->fhv != NULL) fhv_delete_it(bfd->fhv);
  398. if (bfd->lpf != NULL) lpf_delete_it(bfd->lpf);
  399. xfree(bfd);
  400. return;
  401. }
  402. /* eof */