ivecop.c 7.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334
  1. /**************************************************************************
  2. **
  3. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4. **
  5. ** Meschach Library
  6. **
  7. ** This Meschach Library is provided "as is" without any express
  8. ** or implied warranty of any kind with respect to this software.
  9. ** In particular the authors shall not be liable for any direct,
  10. ** indirect, special, incidental or consequential damages arising
  11. ** in any way from use of the software.
  12. **
  13. ** Everyone is granted permission to copy, modify and redistribute this
  14. ** Meschach Library, provided:
  15. ** 1. All copies contain this copyright notice.
  16. ** 2. All modified copies shall carry a notice stating who
  17. ** made the last modification and the date of such modification.
  18. ** 3. No charge is made for this software or works derived from it.
  19. ** This clause shall not be construed as constraining other software
  20. ** distributed on the same medium as this software, nor is a
  21. ** distribution fee considered a charge.
  22. **
  23. ***************************************************************************/
  24. /* ivecop.c */
  25. #include <stdio.h>
  26. #include "matrix.h"
  27. static char rcsid[] = "$Id: ivecop.c,v 1.6 1996/08/20 18:19:21 stewart Exp $";
  28. static char line[MAXLINE];
  29. /* iv_get -- get integer vector -- see also memory.c */
  30. #ifndef ANSI_C
  31. IVEC *iv_get(dim)
  32. int dim;
  33. #else
  34. IVEC *iv_get(int dim)
  35. #endif
  36. {
  37. IVEC *iv;
  38. /* unsigned int i; */
  39. if (dim < 0)
  40. error(E_NEG,"iv_get");
  41. if ((iv=NEW(IVEC)) == IVNULL )
  42. error(E_MEM,"iv_get");
  43. else if (mem_info_is_on()) {
  44. mem_bytes(TYPE_IVEC,0,sizeof(IVEC));
  45. mem_numvar(TYPE_IVEC,1);
  46. }
  47. iv->dim = iv->max_dim = dim;
  48. if ((iv->ive = NEW_A(dim,int)) == (int *)NULL )
  49. error(E_MEM,"iv_get");
  50. else if (mem_info_is_on()) {
  51. mem_bytes(TYPE_IVEC,0,dim*sizeof(int));
  52. }
  53. return (iv);
  54. }
  55. /* iv_free -- returns iv & asoociated memory back to memory heap */
  56. #ifndef ANSI_C
  57. int iv_free(iv)
  58. IVEC *iv;
  59. #else
  60. int iv_free(IVEC *iv)
  61. #endif
  62. {
  63. if ( iv==IVNULL || iv->dim > MAXDIM )
  64. /* don't trust it */
  65. return (-1);
  66. if ( iv->ive == (int *)NULL ) {
  67. if (mem_info_is_on()) {
  68. mem_bytes(TYPE_IVEC,sizeof(IVEC),0);
  69. mem_numvar(TYPE_IVEC,-1);
  70. }
  71. free((char *)iv);
  72. }
  73. else
  74. {
  75. if (mem_info_is_on()) {
  76. mem_bytes(TYPE_IVEC,sizeof(IVEC)+iv->max_dim*sizeof(int),0);
  77. mem_numvar(TYPE_IVEC,-1);
  78. }
  79. free((char *)iv->ive);
  80. free((char *)iv);
  81. }
  82. return (0);
  83. }
  84. /* iv_resize -- returns the IVEC with dimension new_dim
  85. -- iv is set to the zero vector */
  86. #ifndef ANSI_C
  87. IVEC *iv_resize(iv,new_dim)
  88. IVEC *iv;
  89. int new_dim;
  90. #else
  91. IVEC *iv_resize(IVEC *iv, int new_dim)
  92. #endif
  93. {
  94. int i;
  95. if (new_dim < 0)
  96. error(E_NEG,"iv_resize");
  97. if ( ! iv )
  98. return iv_get(new_dim);
  99. if (new_dim == iv->dim)
  100. return iv;
  101. if ( new_dim > iv->max_dim )
  102. {
  103. if (mem_info_is_on()) {
  104. mem_bytes(TYPE_IVEC,iv->max_dim*sizeof(int),
  105. new_dim*sizeof(int));
  106. }
  107. iv->ive = RENEW(iv->ive,new_dim,int);
  108. if ( ! iv->ive )
  109. error(E_MEM,"iv_resize");
  110. iv->max_dim = new_dim;
  111. }
  112. if ( iv->dim <= new_dim )
  113. for ( i = iv->dim; i < new_dim; i++ )
  114. iv->ive[i] = 0;
  115. iv->dim = new_dim;
  116. return iv;
  117. }
  118. /* iv_copy -- copy integer vector in to out
  119. -- out created/resized if necessary */
  120. #ifndef ANSI_C
  121. IVEC *iv_copy(in,out)
  122. IVEC *in, *out;
  123. #else
  124. IVEC *iv_copy(const IVEC *in, IVEC *out)
  125. #endif
  126. {
  127. int i;
  128. if ( ! in )
  129. error(E_NULL,"iv_copy");
  130. out = iv_resize(out,in->dim);
  131. for ( i = 0; i < in->dim; i++ )
  132. out->ive[i] = in->ive[i];
  133. return out;
  134. }
  135. /* iv_move -- move selected pieces of an IVEC
  136. -- moves the length dim0 subvector with initial index i0
  137. to the corresponding subvector of out with initial index i1
  138. -- out is resized if necessary */
  139. #ifndef ANSI_C
  140. IVEC *iv_move(in,i0,dim0,out,i1)
  141. IVEC *in, *out;
  142. int i0, dim0, i1;
  143. #else
  144. IVEC *iv_move(const IVEC *in, int i0, int dim0, IVEC *out, int i1)
  145. #endif
  146. {
  147. if ( ! in )
  148. error(E_NULL,"iv_move");
  149. if ( i0 < 0 || dim0 < 0 || i1 < 0 ||
  150. i0+dim0 > in->dim )
  151. error(E_BOUNDS,"iv_move");
  152. if ( (! out) || i1+dim0 > out->dim )
  153. out = iv_resize(out,i1+dim0);
  154. MEM_COPY(&(in->ive[i0]),&(out->ive[i1]),dim0*sizeof(int));
  155. return out;
  156. }
  157. /* iv_add -- integer vector addition -- may be in-situ */
  158. #ifndef ANSI_C
  159. IVEC *iv_add(iv1,iv2,out)
  160. IVEC *iv1,*iv2,*out;
  161. #else
  162. IVEC *iv_add(const IVEC *iv1, const IVEC *iv2, IVEC *out)
  163. #endif
  164. {
  165. unsigned int i;
  166. int *out_ive, *iv1_ive, *iv2_ive;
  167. if ( iv1==IVNULL || iv2==IVNULL )
  168. error(E_NULL,"iv_add");
  169. if ( iv1->dim != iv2->dim )
  170. error(E_SIZES,"iv_add");
  171. if ( out==IVNULL || out->dim != iv1->dim )
  172. out = iv_resize(out,iv1->dim);
  173. out_ive = out->ive;
  174. iv1_ive = iv1->ive;
  175. iv2_ive = iv2->ive;
  176. for ( i = 0; i < iv1->dim; i++ )
  177. out_ive[i] = iv1_ive[i] + iv2_ive[i];
  178. return (out);
  179. }
  180. /* iv_sub -- integer vector addition -- may be in-situ */
  181. #ifndef ANSI_C
  182. IVEC *iv_sub(iv1,iv2,out)
  183. IVEC *iv1,*iv2,*out;
  184. #else
  185. IVEC *iv_sub(const IVEC *iv1, const IVEC *iv2, IVEC *out)
  186. #endif
  187. {
  188. unsigned int i;
  189. int *out_ive, *iv1_ive, *iv2_ive;
  190. if ( iv1==IVNULL || iv2==IVNULL )
  191. error(E_NULL,"iv_sub");
  192. if ( iv1->dim != iv2->dim )
  193. error(E_SIZES,"iv_sub");
  194. if ( out==IVNULL || out->dim != iv1->dim )
  195. out = iv_resize(out,iv1->dim);
  196. out_ive = out->ive;
  197. iv1_ive = iv1->ive;
  198. iv2_ive = iv2->ive;
  199. for ( i = 0; i < iv1->dim; i++ )
  200. out_ive[i] = iv1_ive[i] - iv2_ive[i];
  201. return (out);
  202. }
  203. #define MAX_STACK 60
  204. /* iv_sort -- sorts vector x, and generates permutation that gives the order
  205. of the components; x = [1.3, 3.7, 0.5] -> [0.5, 1.3, 3.7] and
  206. the permutation is order = [2, 0, 1].
  207. -- if order is NULL on entry then it is ignored
  208. -- the sorted vector x is returned */
  209. #ifndef ANSI_C
  210. IVEC *iv_sort(x, order)
  211. IVEC *x;
  212. PERM *order;
  213. #else
  214. IVEC *iv_sort(IVEC *x, PERM *order)
  215. #endif
  216. {
  217. int *x_ive, tmp, v;
  218. /* int *order_pe; */
  219. int dim, i, j, l, r, tmp_i;
  220. int stack[MAX_STACK], sp;
  221. if ( ! x )
  222. error(E_NULL,"iv_sort");
  223. if ( order != PNULL && order->size != x->dim )
  224. order = px_resize(order, x->dim);
  225. x_ive = x->ive;
  226. dim = x->dim;
  227. if ( order != PNULL )
  228. px_ident(order);
  229. if ( dim <= 1 )
  230. return x;
  231. /* using quicksort algorithm in Sedgewick,
  232. "Algorithms in C", Ch. 9, pp. 118--122 (1990) */
  233. sp = 0;
  234. l = 0; r = dim-1; v = x_ive[0];
  235. for ( ; ; )
  236. {
  237. while ( r > l )
  238. {
  239. /* "i = partition(x_ive,l,r);" */
  240. v = x_ive[r];
  241. i = l-1;
  242. j = r;
  243. for ( ; ; )
  244. {
  245. while ( x_ive[++i] < v )
  246. ;
  247. --j;
  248. while ( x_ive[j] > v && j != 0 )
  249. --j;
  250. if ( i >= j ) break;
  251. tmp = x_ive[i];
  252. x_ive[i] = x_ive[j];
  253. x_ive[j] = tmp;
  254. if ( order != PNULL )
  255. {
  256. tmp_i = order->pe[i];
  257. order->pe[i] = order->pe[j];
  258. order->pe[j] = tmp_i;
  259. }
  260. }
  261. tmp = x_ive[i];
  262. x_ive[i] = x_ive[r];
  263. x_ive[r] = tmp;
  264. if ( order != PNULL )
  265. {
  266. tmp_i = order->pe[i];
  267. order->pe[i] = order->pe[r];
  268. order->pe[r] = tmp_i;
  269. }
  270. if ( i-l > r-i )
  271. { stack[sp++] = l; stack[sp++] = i-1; l = i+1; }
  272. else
  273. { stack[sp++] = i+1; stack[sp++] = r; r = i-1; }
  274. }
  275. /* recursion elimination */
  276. if ( sp == 0 )
  277. break;
  278. r = stack[--sp];
  279. l = stack[--sp];
  280. }
  281. return x;
  282. }