spchfctr.c 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659
  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 Cholesky factorisation code
  26. To be used with sparse.h, sparse.c etc
  27. */
  28. static char rcsid[] = "$Id: spchfctr.c,v 1.5 1996/08/20 19:45:33 stewart Exp $";
  29. #include <stdio.h>
  30. #include <math.h>
  31. #include "sparse2.h"
  32. #ifndef MALLOCDECL
  33. #ifndef ANSI_C
  34. extern char *calloc(), *realloc();
  35. #endif
  36. #endif
  37. /* sprow_ip -- finds the (partial) inner product of a pair of sparse rows
  38. -- uses a "merging" approach & assumes column ordered rows
  39. -- row indices for inner product are all < lim */
  40. #ifndef ANSI_C
  41. static double sprow_ip(row1, row2, lim)
  42. SPROW *row1, *row2;
  43. int lim;
  44. #else
  45. static double sprow_ip(const SPROW *row1, const SPROW *row2, int lim)
  46. #endif
  47. {
  48. int idx1, idx2, len1, len2, tmp;
  49. register row_elt *elts1, *elts2;
  50. register Real sum;
  51. elts1 = row1->elt; elts2 = row2->elt;
  52. len1 = row1->len; len2 = row2->len;
  53. sum = 0.0;
  54. if ( len1 <= 0 || len2 <= 0 )
  55. return 0.0;
  56. if ( elts1->col >= lim || elts2->col >= lim )
  57. return 0.0;
  58. /* use sprow_idx() to speed up inner product where one row is
  59. much longer than the other */
  60. idx1 = idx2 = 0;
  61. if ( len1 > 2*len2 )
  62. {
  63. idx1 = sprow_idx(row1,elts2->col);
  64. idx1 = (idx1 < 0) ? -(idx1+2) : idx1;
  65. if ( idx1 < 0 )
  66. error(E_UNKNOWN,"sprow_ip");
  67. len1 -= idx1;
  68. }
  69. else if ( len2 > 2*len1 )
  70. {
  71. idx2 = sprow_idx(row2,elts1->col);
  72. idx2 = (idx2 < 0) ? -(idx2+2) : idx2;
  73. if ( idx2 < 0 )
  74. error(E_UNKNOWN,"sprow_ip");
  75. len2 -= idx2;
  76. }
  77. if ( len1 <= 0 || len2 <= 0 )
  78. return 0.0;
  79. elts1 = &(elts1[idx1]); elts2 = &(elts2[idx2]);
  80. for ( ; ; ) /* forever do... */
  81. {
  82. if ( (tmp=elts1->col-elts2->col) < 0 )
  83. {
  84. len1--; elts1++;
  85. if ( ! len1 || elts1->col >= lim )
  86. break;
  87. }
  88. else if ( tmp > 0 )
  89. {
  90. len2--; elts2++;
  91. if ( ! len2 || elts2->col >= lim )
  92. break;
  93. }
  94. else
  95. {
  96. sum += elts1->val * elts2->val;
  97. len1--; elts1++;
  98. len2--; elts2++;
  99. if ( ! len1 || ! len2 ||
  100. elts1->col >= lim || elts2->col >= lim )
  101. break;
  102. }
  103. }
  104. return sum;
  105. }
  106. /* sprow_sqr -- returns same as sprow_ip(row, row, lim) */
  107. #ifndef ANSI_C
  108. static double sprow_sqr(row, lim)
  109. SPROW *row;
  110. int lim;
  111. #else
  112. static double sprow_sqr(const SPROW *row, int lim)
  113. #endif
  114. {
  115. register row_elt *elts;
  116. int idx, len;
  117. register Real sum, tmp;
  118. sum = 0.0;
  119. elts = row->elt; len = row->len;
  120. for ( idx = 0; idx < len; idx++, elts++ )
  121. {
  122. if ( elts->col >= lim )
  123. break;
  124. tmp = elts->val;
  125. sum += tmp*tmp;
  126. }
  127. return sum;
  128. }
  129. static int *scan_row = (int *)NULL, *scan_idx = (int *)NULL,
  130. *col_list = (int *)NULL;
  131. static int scan_len = 0;
  132. /* set_scan -- expand scan_row and scan_idx arrays
  133. -- return new length */
  134. #ifndef ANSI_C
  135. int set_scan(new_len)
  136. int new_len;
  137. #else
  138. int set_scan(int new_len)
  139. #endif
  140. {
  141. if ( new_len <= scan_len )
  142. return scan_len;
  143. if ( new_len <= scan_len+5 )
  144. new_len += 5;
  145. /* update scan_len */
  146. scan_len = new_len;
  147. if ( ! scan_row || ! scan_idx || ! col_list )
  148. {
  149. scan_row = (int *)calloc(new_len,sizeof(int));
  150. scan_idx = (int *)calloc(new_len,sizeof(int));
  151. col_list = (int *)calloc(new_len,sizeof(int));
  152. }
  153. else
  154. {
  155. scan_row = (int *)realloc((char *)scan_row,new_len*sizeof(int));
  156. scan_idx = (int *)realloc((char *)scan_idx,new_len*sizeof(int));
  157. col_list = (int *)realloc((char *)col_list,new_len*sizeof(int));
  158. }
  159. if ( ! scan_row || ! scan_idx || ! col_list )
  160. error(E_MEM,"set_scan");
  161. return new_len;
  162. }
  163. /* spCHfactor -- sparse Cholesky factorisation
  164. -- only the lower triangular part of A (incl. diagonal) is used */
  165. #ifndef ANSI_C
  166. SPMAT *spCHfactor(A)
  167. SPMAT *A;
  168. #else
  169. SPMAT *spCHfactor(SPMAT *A)
  170. #endif
  171. {
  172. register int i;
  173. int idx, k, m, minim, n, num_scan, diag_idx, tmp1;
  174. Real pivot, tmp2;
  175. SPROW *r_piv, *r_op;
  176. row_elt *elt_piv, *elt_op, *old_elt;
  177. if ( A == SMNULL )
  178. error(E_NULL,"spCHfactor");
  179. if ( A->m != A->n )
  180. error(E_SQUARE,"spCHfactor");
  181. /* set up access paths if not already done so */
  182. sp_col_access(A);
  183. sp_diag_access(A);
  184. /* printf("spCHfactor() -- checkpoint 1\n"); */
  185. m = A->m; n = A->n;
  186. for ( k = 0; k < m; k++ )
  187. {
  188. r_piv = &(A->row[k]);
  189. if ( r_piv->len > scan_len )
  190. set_scan(r_piv->len);
  191. elt_piv = r_piv->elt;
  192. diag_idx = sprow_idx2(r_piv,k,r_piv->diag);
  193. if ( diag_idx < 0 )
  194. error(E_POSDEF,"spCHfactor");
  195. old_elt = &(elt_piv[diag_idx]);
  196. for ( i = 0; i < r_piv->len; i++ )
  197. {
  198. if ( elt_piv[i].col > k )
  199. break;
  200. col_list[i] = elt_piv[i].col;
  201. scan_row[i] = elt_piv[i].nxt_row;
  202. scan_idx[i] = elt_piv[i].nxt_idx;
  203. }
  204. /* printf("spCHfactor() -- checkpoint 2\n"); */
  205. num_scan = i; /* number of actual entries in scan_row etc. */
  206. /* printf("num_scan = %d\n",num_scan); */
  207. /* set diagonal entry of Cholesky factor */
  208. tmp2 = elt_piv[diag_idx].val - sprow_sqr(r_piv,k);
  209. if ( tmp2 <= 0.0 )
  210. error(E_POSDEF,"spCHfactor");
  211. elt_piv[diag_idx].val = pivot = sqrt(tmp2);
  212. /* now set the k-th column of the Cholesky factors */
  213. /* printf("k = %d\n",k); */
  214. for ( ; ; ) /* forever do... */
  215. {
  216. /* printf("spCHfactor() -- checkpoint 3\n"); */
  217. /* find next row where something (non-trivial) happens
  218. i.e. find min(scan_row) */
  219. /* printf("scan_row: "); */
  220. minim = n;
  221. for ( i = 0; i < num_scan; i++ )
  222. {
  223. tmp1 = scan_row[i];
  224. /* printf("%d ",tmp1); */
  225. minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
  226. }
  227. /* printf("minim = %d\n",minim); */
  228. /* printf("col_list: "); */
  229. /* for ( i = 0; i < num_scan; i++ ) */
  230. /* printf("%d ",col_list[i]); */
  231. /* printf("\n"); */
  232. if ( minim >= n )
  233. break; /* nothing more to do for this column */
  234. r_op = &(A->row[minim]);
  235. elt_op = r_op->elt;
  236. /* set next entry in column k of Cholesky factors */
  237. idx = sprow_idx2(r_op,k,scan_idx[num_scan-1]);
  238. if ( idx < 0 )
  239. { /* fill-in */
  240. sp_set_val(A,minim,k,
  241. -sprow_ip(r_piv,r_op,k)/pivot);
  242. /* in case a realloc() has occurred... */
  243. elt_op = r_op->elt;
  244. /* now set up column access path again */
  245. idx = sprow_idx2(r_op,k,-(idx+2));
  246. tmp1 = old_elt->nxt_row;
  247. old_elt->nxt_row = minim;
  248. r_op->elt[idx].nxt_row = tmp1;
  249. tmp1 = old_elt->nxt_idx;
  250. old_elt->nxt_idx = idx;
  251. r_op->elt[idx].nxt_idx = tmp1;
  252. }
  253. else
  254. elt_op[idx].val = (elt_op[idx].val -
  255. sprow_ip(r_piv,r_op,k))/pivot;
  256. /* printf("spCHfactor() -- checkpoint 4\n"); */
  257. /* remember current element in column k for column chain */
  258. idx = sprow_idx2(r_op,k,idx);
  259. old_elt = &(r_op->elt[idx]);
  260. /* update scan_row */
  261. /* printf("spCHfactor() -- checkpoint 5\n"); */
  262. /* printf("minim = %d\n",minim); */
  263. for ( i = 0; i < num_scan; i++ )
  264. {
  265. if ( scan_row[i] != minim )
  266. continue;
  267. idx = sprow_idx2(r_op,col_list[i],scan_idx[i]);
  268. if ( idx < 0 )
  269. { scan_row[i] = -1; continue; }
  270. scan_row[i] = elt_op[idx].nxt_row;
  271. scan_idx[i] = elt_op[idx].nxt_idx;
  272. /* printf("scan_row[%d] = %d\n",i,scan_row[i]); */
  273. /* printf("scan_idx[%d] = %d\n",i,scan_idx[i]); */
  274. }
  275. }
  276. /* printf("spCHfactor() -- checkpoint 6\n"); */
  277. /* sp_dump(stdout,A); */
  278. /* printf("\n\n\n"); */
  279. }
  280. return A;
  281. }
  282. /* spCHsolve -- solve L.L^T.out=b where L is a sparse matrix,
  283. -- out, b dense vectors
  284. -- returns out; operation may be in-situ */
  285. #ifndef ANSI_C
  286. VEC *spCHsolve(L,b,out)
  287. SPMAT *L;
  288. VEC *b, *out;
  289. #else
  290. VEC *spCHsolve(SPMAT *L, const VEC *b, VEC *out)
  291. #endif
  292. {
  293. int i, j_idx, n, scan_idx, scan_row;
  294. SPROW *row;
  295. row_elt *elt;
  296. Real diag_val, sum, *out_ve;
  297. if ( L == SMNULL || b == VNULL )
  298. error(E_NULL,"spCHsolve");
  299. if ( L->m != L->n )
  300. error(E_SQUARE,"spCHsolve");
  301. if ( b->dim != L->m )
  302. error(E_SIZES,"spCHsolve");
  303. if ( ! L->flag_col )
  304. sp_col_access(L);
  305. if ( ! L->flag_diag )
  306. sp_diag_access(L);
  307. out = v_copy(b,out);
  308. out_ve = out->ve;
  309. /* forward substitution: solve L.x=b for x */
  310. n = L->n;
  311. for ( i = 0; i < n; i++ )
  312. {
  313. sum = out_ve[i];
  314. row = &(L->row[i]);
  315. elt = row->elt;
  316. for ( j_idx = 0; j_idx < row->len; j_idx++, elt++ )
  317. {
  318. if ( elt->col >= i )
  319. break;
  320. sum -= elt->val*out_ve[elt->col];
  321. }
  322. if ( row->diag >= 0 )
  323. out_ve[i] = sum/(row->elt[row->diag].val);
  324. else
  325. error(E_SING,"spCHsolve");
  326. }
  327. /* backward substitution: solve L^T.out = x for out */
  328. for ( i = n-1; i >= 0; i-- )
  329. {
  330. sum = out_ve[i];
  331. row = &(L->row[i]);
  332. /* Note that row->diag >= 0 by above loop */
  333. elt = &(row->elt[row->diag]);
  334. diag_val = elt->val;
  335. /* scan down column */
  336. scan_idx = elt->nxt_idx;
  337. scan_row = elt->nxt_row;
  338. while ( scan_row >= 0 /* && scan_idx >= 0 */ )
  339. {
  340. row = &(L->row[scan_row]);
  341. elt = &(row->elt[scan_idx]);
  342. sum -= elt->val*out_ve[scan_row];
  343. scan_idx = elt->nxt_idx;
  344. scan_row = elt->nxt_row;
  345. }
  346. out_ve[i] = sum/diag_val;
  347. }
  348. return out;
  349. }
  350. /* spICHfactor -- sparse Incomplete Cholesky factorisation
  351. -- does a Cholesky factorisation assuming NO FILL-IN
  352. -- as for spCHfactor(), only the lower triangular part of A is used */
  353. #ifndef ANSI_C
  354. SPMAT *spICHfactor(A)
  355. SPMAT *A;
  356. #else
  357. SPMAT *spICHfactor(SPMAT *A)
  358. #endif
  359. {
  360. int k, m, n, nxt_row, nxt_idx, diag_idx;
  361. Real pivot, tmp2;
  362. SPROW *r_piv, *r_op;
  363. row_elt *elt_piv, *elt_op;
  364. if ( A == SMNULL )
  365. error(E_NULL,"spICHfactor");
  366. if ( A->m != A->n )
  367. error(E_SQUARE,"spICHfactor");
  368. /* set up access paths if not already done so */
  369. if ( ! A->flag_col )
  370. sp_col_access(A);
  371. if ( ! A->flag_diag )
  372. sp_diag_access(A);
  373. m = A->m; n = A->n;
  374. for ( k = 0; k < m; k++ )
  375. {
  376. r_piv = &(A->row[k]);
  377. diag_idx = r_piv->diag;
  378. if ( diag_idx < 0 )
  379. error(E_POSDEF,"spICHfactor");
  380. elt_piv = r_piv->elt;
  381. /* set diagonal entry of Cholesky factor */
  382. tmp2 = elt_piv[diag_idx].val - sprow_sqr(r_piv,k);
  383. if ( tmp2 <= 0.0 )
  384. error(E_POSDEF,"spICHfactor");
  385. elt_piv[diag_idx].val = pivot = sqrt(tmp2);
  386. /* find next row where something (non-trivial) happens */
  387. nxt_row = elt_piv[diag_idx].nxt_row;
  388. nxt_idx = elt_piv[diag_idx].nxt_idx;
  389. /* now set the k-th column of the Cholesky factors */
  390. while ( nxt_row >= 0 && nxt_idx >= 0 )
  391. {
  392. /* nxt_row and nxt_idx give next next row (& index)
  393. of the entry to be modified */
  394. r_op = &(A->row[nxt_row]);
  395. elt_op = r_op->elt;
  396. elt_op[nxt_idx].val = (elt_op[nxt_idx].val -
  397. sprow_ip(r_piv,r_op,k))/pivot;
  398. nxt_row = elt_op[nxt_idx].nxt_row;
  399. nxt_idx = elt_op[nxt_idx].nxt_idx;
  400. }
  401. }
  402. return A;
  403. }
  404. /* spCHsymb -- symbolic sparse Cholesky factorisation
  405. -- does NOT do any floating point arithmetic; just sets up the structure
  406. -- only the lower triangular part of A (incl. diagonal) is used */
  407. #ifndef ANSI_C
  408. SPMAT *spCHsymb(A)
  409. SPMAT *A;
  410. #else
  411. SPMAT *spCHsymb(SPMAT *A)
  412. #endif
  413. {
  414. register int i;
  415. int idx, k, m, minim, n, num_scan, diag_idx, tmp1;
  416. SPROW *r_piv, *r_op;
  417. row_elt *elt_piv, *elt_op, *old_elt;
  418. if ( A == SMNULL )
  419. error(E_NULL,"spCHsymb");
  420. if ( A->m != A->n )
  421. error(E_SQUARE,"spCHsymb");
  422. /* set up access paths if not already done so */
  423. if ( ! A->flag_col )
  424. sp_col_access(A);
  425. if ( ! A->flag_diag )
  426. sp_diag_access(A);
  427. /* printf("spCHsymb() -- checkpoint 1\n"); */
  428. m = A->m; n = A->n;
  429. for ( k = 0; k < m; k++ )
  430. {
  431. r_piv = &(A->row[k]);
  432. if ( r_piv->len > scan_len )
  433. set_scan(r_piv->len);
  434. elt_piv = r_piv->elt;
  435. diag_idx = sprow_idx2(r_piv,k,r_piv->diag);
  436. if ( diag_idx < 0 )
  437. error(E_POSDEF,"spCHsymb");
  438. old_elt = &(elt_piv[diag_idx]);
  439. for ( i = 0; i < r_piv->len; i++ )
  440. {
  441. if ( elt_piv[i].col > k )
  442. break;
  443. col_list[i] = elt_piv[i].col;
  444. scan_row[i] = elt_piv[i].nxt_row;
  445. scan_idx[i] = elt_piv[i].nxt_idx;
  446. }
  447. /* printf("spCHsymb() -- checkpoint 2\n"); */
  448. num_scan = i; /* number of actual entries in scan_row etc. */
  449. /* printf("num_scan = %d\n",num_scan); */
  450. /* now set the k-th column of the Cholesky factors */
  451. /* printf("k = %d\n",k); */
  452. for ( ; ; ) /* forever do... */
  453. {
  454. /* printf("spCHsymb() -- checkpoint 3\n"); */
  455. /* find next row where something (non-trivial) happens
  456. i.e. find min(scan_row) */
  457. minim = n;
  458. for ( i = 0; i < num_scan; i++ )
  459. {
  460. tmp1 = scan_row[i];
  461. /* printf("%d ",tmp1); */
  462. minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
  463. }
  464. if ( minim >= n )
  465. break; /* nothing more to do for this column */
  466. r_op = &(A->row[minim]);
  467. elt_op = r_op->elt;
  468. /* set next entry in column k of Cholesky factors */
  469. idx = sprow_idx2(r_op,k,scan_idx[num_scan-1]);
  470. if ( idx < 0 )
  471. { /* fill-in */
  472. sp_set_val(A,minim,k,0.0);
  473. /* in case a realloc() has occurred... */
  474. elt_op = r_op->elt;
  475. /* now set up column access path again */
  476. idx = sprow_idx2(r_op,k,-(idx+2));
  477. tmp1 = old_elt->nxt_row;
  478. old_elt->nxt_row = minim;
  479. r_op->elt[idx].nxt_row = tmp1;
  480. tmp1 = old_elt->nxt_idx;
  481. old_elt->nxt_idx = idx;
  482. r_op->elt[idx].nxt_idx = tmp1;
  483. }
  484. /* printf("spCHsymb() -- checkpoint 4\n"); */
  485. /* remember current element in column k for column chain */
  486. idx = sprow_idx2(r_op,k,idx);
  487. old_elt = &(r_op->elt[idx]);
  488. /* update scan_row */
  489. /* printf("spCHsymb() -- checkpoint 5\n"); */
  490. /* printf("minim = %d\n",minim); */
  491. for ( i = 0; i < num_scan; i++ )
  492. {
  493. if ( scan_row[i] != minim )
  494. continue;
  495. idx = sprow_idx2(r_op,col_list[i],scan_idx[i]);
  496. if ( idx < 0 )
  497. { scan_row[i] = -1; continue; }
  498. scan_row[i] = elt_op[idx].nxt_row;
  499. scan_idx[i] = elt_op[idx].nxt_idx;
  500. /* printf("scan_row[%d] = %d\n",i,scan_row[i]); */
  501. /* printf("scan_idx[%d] = %d\n",i,scan_idx[i]); */
  502. }
  503. }
  504. /* printf("spCHsymb() -- checkpoint 6\n"); */
  505. }
  506. return A;
  507. }
  508. /* comp_AAT -- compute A.A^T where A is a given sparse matrix */
  509. #ifndef ANSI_C
  510. SPMAT *comp_AAT(A)
  511. SPMAT *A;
  512. #else
  513. SPMAT *comp_AAT(SPMAT *A)
  514. #endif
  515. {
  516. SPMAT *AAT;
  517. SPROW *r, *r2;
  518. row_elt *elts, *elts2;
  519. int i, idx, idx2, j, m, minim, n, num_scan, tmp1;
  520. Real ip;
  521. if ( ! A )
  522. error(E_NULL,"comp_AAT");
  523. m = A->m; n = A->n;
  524. /* set up column access paths */
  525. if ( ! A->flag_col )
  526. sp_col_access(A);
  527. AAT = sp_get(m,m,10);
  528. for ( i = 0; i < m; i++ )
  529. {
  530. /* initialisation */
  531. r = &(A->row[i]);
  532. elts = r->elt;
  533. /* set up scan lists for this row */
  534. if ( r->len > scan_len )
  535. set_scan(r->len);
  536. for ( j = 0; j < r->len; j++ )
  537. {
  538. col_list[j] = elts[j].col;
  539. scan_row[j] = elts[j].nxt_row;
  540. scan_idx[j] = elts[j].nxt_idx;
  541. }
  542. num_scan = r->len;
  543. /* scan down the rows for next non-zero not
  544. associated with a diagonal entry */
  545. for ( ; ; )
  546. {
  547. minim = m;
  548. for ( idx = 0; idx < num_scan; idx++ )
  549. {
  550. tmp1 = scan_row[idx];
  551. minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
  552. }
  553. if ( minim >= m )
  554. break;
  555. r2 = &(A->row[minim]);
  556. if ( minim > i )
  557. {
  558. ip = sprow_ip(r,r2,n);
  559. sp_set_val(AAT,minim,i,ip);
  560. sp_set_val(AAT,i,minim,ip);
  561. }
  562. /* update scan entries */
  563. elts2 = r2->elt;
  564. for ( idx = 0; idx < num_scan; idx++ )
  565. {
  566. if ( scan_row[idx] != minim || scan_idx[idx] < 0 )
  567. continue;
  568. idx2 = scan_idx[idx];
  569. scan_row[idx] = elts2[idx2].nxt_row;
  570. scan_idx[idx] = elts2[idx2].nxt_idx;
  571. }
  572. }
  573. /* set the diagonal entry */
  574. sp_set_val(AAT,i,i,sprow_sqr(r,n));
  575. }
  576. return AAT;
  577. }