spswap.c 7.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326
  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. /*
  25. Sparse matrix swap and permutation routines
  26. Modified Mon 09th Nov 1992, 08:50:54 PM
  27. to use Karen George's suggestion to use unordered rows
  28. */
  29. static char rcsid[] = "$Id: spswap.c,v 1.3 1994/01/13 05:44:43 des Exp $";
  30. #include <stdio.h>
  31. #include <math.h>
  32. #include "sparse2.h"
  33. #define btos(x) ((x) ? "TRUE" : "FALSE")
  34. /* scan_to -- updates scan (int) vectors to point to the last row in each
  35. column with row # <= max_row, if any */
  36. #ifndef ANSI_C
  37. void scan_to(A, scan_row, scan_idx, col_list, max_row)
  38. SPMAT *A;
  39. IVEC *scan_row, *scan_idx, *col_list;
  40. int max_row;
  41. #else
  42. void scan_to(SPMAT *A, IVEC *scan_row, IVEC *scan_idx, IVEC *col_list,
  43. int max_row)
  44. #endif
  45. {
  46. int col, idx, j_idx, row_num;
  47. SPROW *r;
  48. row_elt *e;
  49. if ( ! A || ! scan_row || ! scan_idx || ! col_list )
  50. error(E_NULL,"scan_to");
  51. if ( scan_row->dim != scan_idx->dim || scan_idx->dim != col_list->dim )
  52. error(E_SIZES,"scan_to");
  53. if ( max_row < 0 )
  54. return;
  55. if ( ! A->flag_col )
  56. sp_col_access(A);
  57. for ( j_idx = 0; j_idx < scan_row->dim; j_idx++ )
  58. {
  59. row_num = scan_row->ive[j_idx];
  60. idx = scan_idx->ive[j_idx];
  61. col = col_list->ive[j_idx];
  62. if ( col < 0 || col >= A->n )
  63. error(E_BOUNDS,"scan_to");
  64. if ( row_num < 0 )
  65. {
  66. idx = col;
  67. continue;
  68. }
  69. r = &(A->row[row_num]);
  70. if ( idx < 0 )
  71. error(E_INTERN,"scan_to");
  72. e = &(r->elt[idx]);
  73. if ( e->col != col )
  74. error(E_INTERN,"scan_to");
  75. if ( idx < 0 )
  76. {
  77. printf("scan_to: row_num = %d, idx = %d, col = %d\n",
  78. row_num, idx, col);
  79. error(E_INTERN,"scan_to");
  80. }
  81. /* if ( e->nxt_row <= max_row )
  82. chase_col(A, col, &row_num, &idx, max_row); */
  83. while ( e->nxt_row >= 0 && e->nxt_row <= max_row )
  84. {
  85. row_num = e->nxt_row;
  86. idx = e->nxt_idx;
  87. e = &(A->row[row_num].elt[idx]);
  88. }
  89. /* printf("scan_to: computed j_idx = %d, row_num = %d, idx = %d\n",
  90. j_idx, row_num, idx); */
  91. scan_row->ive[j_idx] = row_num;
  92. scan_idx->ive[j_idx] = idx;
  93. }
  94. }
  95. /* patch_col -- patches column access paths for fill-in */
  96. #ifndef ANSI_C
  97. void patch_col(A, col, old_row, old_idx, row_num, idx)
  98. SPMAT *A;
  99. int col, old_row, old_idx, row_num, idx;
  100. #else
  101. void patch_col(SPMAT *A, int col, int old_row, int old_idx, int row_num,
  102. int idx)
  103. #endif
  104. {
  105. SPROW *r;
  106. row_elt *e;
  107. if ( old_row >= 0 )
  108. {
  109. r = &(A->row[old_row]);
  110. old_idx = sprow_idx2(r,col,old_idx);
  111. e = &(r->elt[old_idx]);
  112. e->nxt_row = row_num;
  113. e->nxt_idx = idx;
  114. }
  115. else
  116. {
  117. A->start_row[col] = row_num;
  118. A->start_idx[col] = idx;
  119. }
  120. }
  121. /* chase_col -- chases column access path in column col, starting with
  122. row_num and idx, to find last row # in this column <= max_row
  123. -- row_num is returned; idx is also set by this routine
  124. -- assumes that the column access paths (possibly without the
  125. nxt_idx fields) are set up */
  126. #ifndef ANSI_C
  127. row_elt *chase_col(A, col, row_num, idx, max_row)
  128. SPMAT *A;
  129. int col, *row_num, *idx, max_row;
  130. #else
  131. row_elt *chase_col(const SPMAT *A, int col, int *row_num, int *idx,
  132. int max_row)
  133. #endif
  134. {
  135. int old_idx, old_row, tmp_idx, tmp_row;
  136. SPROW *r;
  137. row_elt *e;
  138. if ( col < 0 || col >= A->n )
  139. error(E_BOUNDS,"chase_col");
  140. tmp_row = *row_num;
  141. if ( tmp_row < 0 )
  142. {
  143. if ( A->start_row[col] > max_row )
  144. {
  145. tmp_row = -1;
  146. tmp_idx = col;
  147. return (row_elt *)NULL;
  148. }
  149. else
  150. {
  151. tmp_row = A->start_row[col];
  152. tmp_idx = A->start_idx[col];
  153. }
  154. }
  155. else
  156. tmp_idx = *idx;
  157. old_row = tmp_row;
  158. old_idx = tmp_idx;
  159. while ( tmp_row >= 0 && tmp_row < max_row )
  160. {
  161. r = &(A->row[tmp_row]);
  162. /* tmp_idx = sprow_idx2(r,col,tmp_idx); */
  163. if ( tmp_idx < 0 || tmp_idx >= r->len ||
  164. r->elt[tmp_idx].col != col )
  165. {
  166. #ifdef DEBUG
  167. printf("chase_col:error: col = %d, row # = %d, idx = %d\n",
  168. col, tmp_row, tmp_idx);
  169. printf("chase_col:error: old_row = %d, old_idx = %d\n",
  170. old_row, old_idx);
  171. printf("chase_col:error: A =\n");
  172. sp_dump(stdout,A);
  173. #endif
  174. error(E_INTERN,"chase_col");
  175. }
  176. e = &(r->elt[tmp_idx]);
  177. old_row = tmp_row;
  178. old_idx = tmp_idx;
  179. tmp_row = e->nxt_row;
  180. tmp_idx = e->nxt_idx;
  181. }
  182. if ( old_row > max_row )
  183. {
  184. old_row = -1;
  185. old_idx = col;
  186. e = (row_elt *)NULL;
  187. }
  188. else if ( tmp_row <= max_row && tmp_row >= 0 )
  189. {
  190. old_row = tmp_row;
  191. old_idx = tmp_idx;
  192. }
  193. *row_num = old_row;
  194. if ( old_row >= 0 )
  195. *idx = old_idx;
  196. else
  197. *idx = col;
  198. return e;
  199. }
  200. /* chase_past -- as for chase_col except that we want the first
  201. row whose row # >= min_row; -1 indicates no such row */
  202. #ifndef ANSI_C
  203. row_elt *chase_past(A, col, row_num, idx, min_row)
  204. SPMAT *A;
  205. int col, *row_num, *idx, min_row;
  206. #else
  207. row_elt *chase_past(const SPMAT *A, int col, int *row_num, int *idx,
  208. int min_row)
  209. #endif
  210. {
  211. SPROW *r;
  212. row_elt *e;
  213. int tmp_idx, tmp_row;
  214. tmp_row = *row_num;
  215. tmp_idx = *idx;
  216. chase_col(A,col,&tmp_row,&tmp_idx,min_row);
  217. if ( tmp_row < 0 ) /* use A->start_row[..] etc. */
  218. {
  219. if ( A->start_row[col] < 0 )
  220. tmp_row = -1;
  221. else
  222. {
  223. tmp_row = A->start_row[col];
  224. tmp_idx = A->start_idx[col];
  225. }
  226. }
  227. else if ( tmp_row < min_row )
  228. {
  229. r = &(A->row[tmp_row]);
  230. if ( tmp_idx < 0 || tmp_idx >= r->len ||
  231. r->elt[tmp_idx].col != col )
  232. error(E_INTERN,"chase_past");
  233. tmp_row = r->elt[tmp_idx].nxt_row;
  234. tmp_idx = r->elt[tmp_idx].nxt_idx;
  235. }
  236. *row_num = tmp_row;
  237. *idx = tmp_idx;
  238. if ( tmp_row < 0 )
  239. e = (row_elt *)NULL;
  240. else
  241. {
  242. if ( tmp_idx < 0 || tmp_idx >= A->row[tmp_row].len ||
  243. A->row[tmp_row].elt[tmp_idx].col != col )
  244. error(E_INTERN,"bump_col");
  245. e = &(A->row[tmp_row].elt[tmp_idx]);
  246. }
  247. return e;
  248. }
  249. /* bump_col -- move along to next nonzero entry in column col after row_num
  250. -- update row_num and idx */
  251. #ifndef ANSI_C
  252. row_elt *bump_col(A, col, row_num, idx)
  253. SPMAT *A;
  254. int col, *row_num, *idx;
  255. #else
  256. row_elt *bump_col(const SPMAT *A, int col, int *row_num, int *idx)
  257. #endif
  258. {
  259. SPROW *r;
  260. row_elt *e;
  261. int tmp_row, tmp_idx;
  262. tmp_row = *row_num;
  263. tmp_idx = *idx;
  264. /* printf("bump_col: col = %d, row# = %d, idx = %d\n",
  265. col, *row_num, *idx); */
  266. if ( tmp_row < 0 )
  267. {
  268. tmp_row = A->start_row[col];
  269. tmp_idx = A->start_idx[col];
  270. }
  271. else
  272. {
  273. r = &(A->row[tmp_row]);
  274. if ( tmp_idx < 0 || tmp_idx >= r->len ||
  275. r->elt[tmp_idx].col != col )
  276. error(E_INTERN,"bump_col");
  277. e = &(r->elt[tmp_idx]);
  278. tmp_row = e->nxt_row;
  279. tmp_idx = e->nxt_idx;
  280. }
  281. if ( tmp_row < 0 )
  282. {
  283. e = (row_elt *)NULL;
  284. tmp_idx = col;
  285. }
  286. else
  287. {
  288. if ( tmp_idx < 0 || tmp_idx >= A->row[tmp_row].len ||
  289. A->row[tmp_row].elt[tmp_idx].col != col )
  290. error(E_INTERN,"bump_col");
  291. e = &(A->row[tmp_row].elt[tmp_idx]);
  292. }
  293. *row_num = tmp_row;
  294. *idx = tmp_idx;
  295. return e;
  296. }