Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/colamd.cc @ 7924:4976f66d469b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
miscellaneous cleanup
author John W. Eaton <jwe@octave.org>
date Fri, 11 Jul 2008 17:59:28 -0400
parents b166043585a8
children 25bc2d31e1bf
comparison
equal deleted inserted replaced
7923:c3d21b9b94b6 7924:4976f66d469b
51 #define SYMAMD_NAME(name) symamd ## name 51 #define SYMAMD_NAME(name) symamd ## name
52 #endif 52 #endif
53 53
54 // The symmetric column elimination tree code take from the Davis LDL code. 54 // The symmetric column elimination tree code take from the Davis LDL code.
55 // Copyright given elsewhere in this file. 55 // Copyright given elsewhere in this file.
56 static 56 static void
57 void symetree (const octave_idx_type *ridx, const octave_idx_type *cidx, 57 symetree (const octave_idx_type *ridx, const octave_idx_type *cidx,
58 octave_idx_type *Parent, octave_idx_type *P, octave_idx_type n) 58 octave_idx_type *Parent, octave_idx_type *P, octave_idx_type n)
59 { 59 {
60 OCTAVE_LOCAL_BUFFER (octave_idx_type, Flag, n); 60 OCTAVE_LOCAL_BUFFER (octave_idx_type, Flag, n);
61 OCTAVE_LOCAL_BUFFER (octave_idx_type, Pinv, (P ? n : 0)); 61 OCTAVE_LOCAL_BUFFER (octave_idx_type, Pinv, (P ? n : 0));
62 if (P) 62 if (P)
63 // If P is present then compute Pinv, the inverse of P 63 // If P is present then compute Pinv, the inverse of P
89 } 89 }
90 } 90 }
91 } 91 }
92 92
93 // The elimination tree post-ordering code below is taken from SuperLU 93 // The elimination tree post-ordering code below is taken from SuperLU
94 static inline 94 static inline octave_idx_type
95 octave_idx_type make_set (octave_idx_type i, octave_idx_type *pp) 95 make_set (octave_idx_type i, octave_idx_type *pp)
96 { 96 {
97 pp[i] = i; 97 pp[i] = i;
98 return i; 98 return i;
99 } 99 }
100 100
101 static inline 101 static inline octave_idx_type
102 octave_idx_type link (octave_idx_type s, octave_idx_type t, octave_idx_type *pp) 102 link (octave_idx_type s, octave_idx_type t, octave_idx_type *pp)
103 { 103 {
104 pp[s] = t; 104 pp[s] = t;
105 return t; 105 return t;
106 } 106 }
107 107
108 static inline 108 static inline octave_idx_type
109 octave_idx_type find (octave_idx_type i, octave_idx_type *pp) 109 find (octave_idx_type i, octave_idx_type *pp)
110 { 110 {
111 register octave_idx_type p, gp; 111 register octave_idx_type p, gp;
112 112
113 p = pp[i]; 113 p = pp[i];
114 gp = pp[p]; 114 gp = pp[p];
115 while (gp != p) { 115
116 pp[i] = gp; 116 while (gp != p)
117 i = gp; 117 {
118 p = pp[i]; 118 pp[i] = gp;
119 gp = pp[p]; 119 i = gp;
120 } 120 p = pp[i];
121 return (p); 121 gp = pp[p];
122 } 122 }
123 123
124 static 124 return p;
125 octave_idx_type etdfs (octave_idx_type v, octave_idx_type *first_kid, 125 }
126 octave_idx_type *next_kid, octave_idx_type *post, 126
127 octave_idx_type postnum) 127 static octave_idx_type
128 { 128 etdfs (octave_idx_type v, octave_idx_type *first_kid,
129 for (octave_idx_type w = first_kid[v]; w != -1; w = next_kid[w]) { 129 octave_idx_type *next_kid, octave_idx_type *post,
130 octave_idx_type postnum)
131 {
132 for (octave_idx_type w = first_kid[v]; w != -1; w = next_kid[w])
130 postnum = etdfs (w, first_kid, next_kid, post, postnum); 133 postnum = etdfs (w, first_kid, next_kid, post, postnum);
131 } 134
132 post[postnum++] = v; 135 post[postnum++] = v;
133 136
134 return postnum; 137 return postnum;
135 } 138 }
136 139
137 static 140 static void
138 void TreePostorder(octave_idx_type n, octave_idx_type *parent, octave_idx_type *post) 141 tree_postorder (octave_idx_type n, octave_idx_type *parent,
142 octave_idx_type *post)
139 { 143 {
140 // Allocate storage for working arrays and results 144 // Allocate storage for working arrays and results
141 OCTAVE_LOCAL_BUFFER (octave_idx_type, first_kid, n+1); 145 OCTAVE_LOCAL_BUFFER (octave_idx_type, first_kid, n+1);
142 OCTAVE_LOCAL_BUFFER (octave_idx_type, next_kid, n+1); 146 OCTAVE_LOCAL_BUFFER (octave_idx_type, next_kid, n+1);
143 147
144 // Set up structure describing children 148 // Set up structure describing children
145 for (octave_idx_type v = 0; v <= n; first_kid[v++] = -1); 149 for (octave_idx_type v = 0; v <= n; first_kid[v++] = -1)
150 /* do nothing */;
151
146 for (octave_idx_type v = n-1; v >= 0; v--) 152 for (octave_idx_type v = n-1; v >= 0; v--)
147 { 153 {
148 octave_idx_type dad = parent[v]; 154 octave_idx_type dad = parent[v];
149 next_kid[v] = first_kid[dad]; 155 next_kid[v] = first_kid[dad];
150 first_kid[dad] = v; 156 first_kid[dad] = v;
152 158
153 // Depth-first search from dummy root vertex #n 159 // Depth-first search from dummy root vertex #n
154 etdfs (n, first_kid, next_kid, post, 0); 160 etdfs (n, first_kid, next_kid, post, 0);
155 } 161 }
156 162
157 static 163 static void
158 void coletree (const octave_idx_type *ridx, const octave_idx_type *colbeg, 164 coletree (const octave_idx_type *ridx, const octave_idx_type *colbeg,
159 octave_idx_type *colend, octave_idx_type *parent, 165 octave_idx_type *colend, octave_idx_type *parent,
160 octave_idx_type nr, octave_idx_type nc) 166 octave_idx_type nr, octave_idx_type nc)
161 { 167 {
162 OCTAVE_LOCAL_BUFFER (octave_idx_type, root, nc); 168 OCTAVE_LOCAL_BUFFER (octave_idx_type, root, nc);
163 OCTAVE_LOCAL_BUFFER (octave_idx_type, pp, nc); 169 OCTAVE_LOCAL_BUFFER (octave_idx_type, pp, nc);
164 OCTAVE_LOCAL_BUFFER (octave_idx_type, firstcol, nr); 170 OCTAVE_LOCAL_BUFFER (octave_idx_type, firstcol, nr);
165 171
166 // Compute firstcol[row] = first nonzero column in row 172 // Compute firstcol[row] = first nonzero column in row
167 for (octave_idx_type row = 0; row < nr; firstcol[row++] = nc); 173 for (octave_idx_type row = 0; row < nr; firstcol[row++] = nc)
174 /* do nothing */;
175
168 for (octave_idx_type col = 0; col < nc; col++) 176 for (octave_idx_type col = 0; col < nc; col++)
169 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++) 177 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++)
170 { 178 {
171 octave_idx_type row = ridx[p]; 179 octave_idx_type row = ridx[p];
172 if (firstcol[row] > col) 180 if (firstcol[row] > col)
398 } 406 }
399 407
400 coletree (ridx, colbeg, colend, etree, n_row, n_col); 408 coletree (ridx, colbeg, colend, etree, n_row, n_col);
401 409
402 // Calculate the tree post-ordering 410 // Calculate the tree post-ordering
403 TreePostorder (n_col, etree, colbeg); 411 tree_postorder (n_col, etree, colbeg);
404 412
405 // return the permutation vector 413 // return the permutation vector
406 NDArray out_perm (dim_vector (1, n_col)); 414 NDArray out_perm (dim_vector (1, n_col));
407 for (octave_idx_type i = 0; i < n_col; i++) 415 for (octave_idx_type i = 0; i < n_col; i++)
408 out_perm(i) = p [colbeg [i]] + 1; 416 out_perm(i) = p [colbeg [i]] + 1;
409 417
410 retval (0) = out_perm; 418 retval(0) = out_perm;
411 419
412 // print stats if spumoni > 0 420 // print stats if spumoni > 0
413 if (spumoni > 0) 421 if (spumoni > 0)
414 COLAMD_NAME (_report) (stats) ; 422 COLAMD_NAME (_report) (stats) ;
415 423
595 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1); 603 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
596 symetree (ridx, cidx, etree, perm, n_col); 604 symetree (ridx, cidx, etree, perm, n_col);
597 605
598 // Calculate the tree post-ordering 606 // Calculate the tree post-ordering
599 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1); 607 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
600 TreePostorder (n_col, etree, post); 608 tree_postorder (n_col, etree, post);
601 609
602 // return the permutation vector 610 // return the permutation vector
603 NDArray out_perm (dim_vector (1, n_col)); 611 NDArray out_perm (dim_vector (1, n_col));
604 for (octave_idx_type i = 0; i < n_col; i++) 612 for (octave_idx_type i = 0; i < n_col; i++)
605 out_perm(i) = perm [post [i]] + 1; 613 out_perm(i) = perm [post [i]] + 1;
606 614
607 retval (0) = out_perm; 615 retval(0) = out_perm;
608 616
609 // print stats if spumoni > 0 617 // print stats if spumoni > 0
610 if (spumoni > 0) 618 if (spumoni > 0)
611 SYMAMD_NAME (_report) (stats) ; 619 SYMAMD_NAME (_report) (stats) ;
612 620
705 { 713 {
706 error ("etree: second argument must be a string"); 714 error ("etree: second argument must be a string");
707 return retval; 715 return retval;
708 } 716 }
709 } 717 }
718
710 // column elimination tree post-ordering (reuse variables) 719 // column elimination tree post-ordering (reuse variables)
711 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1); 720 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
712
713 721
714 if (is_sym) 722 if (is_sym)
715 { 723 {
716 if (n_row != n_col) 724 if (n_row != n_col)
717 { 725 {
718 error ("etree: matrix is marked as symmetric, but not square"); 726 error ("etree: matrix is marked as symmetric, but not square");
719 return retval; 727 return retval;
720 } 728 }
729
721 symetree (ridx, cidx, etree, 0, n_col); 730 symetree (ridx, cidx, etree, 0, n_col);
722 } 731 }
723 else 732 else
724 { 733 {
725 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col); 734 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col);
737 NDArray tree (dim_vector (1, n_col)); 746 NDArray tree (dim_vector (1, n_col));
738 for (octave_idx_type i = 0; i < n_col; i++) 747 for (octave_idx_type i = 0; i < n_col; i++)
739 // We flag a root with n_col while Matlab does it with zero 748 // We flag a root with n_col while Matlab does it with zero
740 // Convert for matlab compatiable output 749 // Convert for matlab compatiable output
741 if (etree[i] == n_col) 750 if (etree[i] == n_col)
742 tree (i) = 0; 751 tree(i) = 0;
743 else 752 else
744 tree (i) = etree[i] + 1; 753 tree(i) = etree[i] + 1;
745 754
746 retval (0) = tree; 755 retval(0) = tree;
747 756
748 if (nargout == 2) 757 if (nargout == 2)
749 { 758 {
750 // Calculate the tree post-ordering 759 // Calculate the tree post-ordering
751 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1); 760 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
752 TreePostorder (n_col, etree, post); 761 tree_postorder (n_col, etree, post);
753 762
754 NDArray postorder (dim_vector (1, n_col)); 763 NDArray postorder (dim_vector (1, n_col));
755 for (octave_idx_type i = 0; i < n_col; i++) 764 for (octave_idx_type i = 0; i < n_col; i++)
756 postorder (i) = post[i] + 1; 765 postorder(i) = post[i] + 1;
757 766
758 retval (1) = postorder; 767 retval(1) = postorder;
759 } 768 }
760 } 769 }
761 770
762 return retval; 771 return retval;
763 } 772 }