amd_preprocess source code

1 /* ========================================================================= */
2 /* === AMD_preprocess ====================================================== */
3 /* ========================================================================= */
4 
5 /* ------------------------------------------------------------------------- */
6 /* AMD, Copyright (c) Timothy A. Davis,					     */
7 /* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
8 /* email: DrTimothyAldenDavis@gmail.com                                      */
9 /* ------------------------------------------------------------------------- */
10 
11 /* Sorts, removes duplicate entries, and transposes from the nonzero pattern of
12  * a column-form matrix A, to obtain the matrix R.  The input matrix can have
13  * duplicate entries and/or unsorted columns (AMD_valid (n,Ap,Ai) must not be
14  * AMD_INVALID).
15  *
16  * This input condition is NOT checked.  This routine is not user-callable.
17  */
18 
19 module amd_preprocess;
20 
21 nothrow @nogc extern(C):
22 
23 import amd_internal;
24 import amd_valid;
25 import amd;
26 
27 /* ========================================================================= */
28 /* === AMD_preprocess ====================================================== */
29 /* ========================================================================= */
30 
31 /* AMD_preprocess does not check its input for errors or allocate workspace.
32  * On input, the condition (AMD_valid (n,n,Ap,Ai) != AMD_INVALID) must hold.
33  */
34 
35 void AMD_preprocess
36 (
37     Int n,		/* input matrix: A is n-by-n */
38     const Int * Ap,	/* size n+1 */
39     const Int * Ai,	/* size nz = Ap [n] */
40 
41     /* output matrix R: */
42     Int * Rp,		/* size n+1 */
43     Int * Ri,		/* size nz (or less, if duplicates present) */
44 
45     Int * W,		/* workspace of size n */
46     Int * Flag	/* workspace of size n */
47 )
48 {
49 
50     /* --------------------------------------------------------------------- */
51     /* local variables */
52     /* --------------------------------------------------------------------- */
53 
54     Int i, j, p, p2 ;
55 
56     ASSERT (AMD_valid (n, n, Ap, Ai) != AMD_INVALID) ;
57 
58     /* --------------------------------------------------------------------- */
59     /* count the entries in each row of A (excluding duplicates) */
60     /* --------------------------------------------------------------------- */
61 
62     for (i = 0 ; i < n ; i++)
63     {
64 	W [i] = 0 ;		/* # of nonzeros in row i (excl duplicates) */
65 	Flag [i] = EMPTY ;	/* Flag [i] = j if i appears in column j */
66     }
67     for (j = 0 ; j < n ; j++)
68     {
69 	p2 = Ap [j+1] ;
70 	for (p = Ap [j] ; p < p2 ; p++)
71 	{
72 	    i = Ai [p] ;
73 	    if (Flag [i] != j)
74 	    {
75 		/* row index i has not yet appeared in column j */
76 		W [i]++ ;	    /* one more entry in row i */
77 		Flag [i] = j ;	    /* flag row index i as appearing in col j*/
78 	    }
79 	}
80     }
81 
82     /* --------------------------------------------------------------------- */
83     /* compute the row pointers for R */
84     /* --------------------------------------------------------------------- */
85 
86     Rp [0] = 0 ;
87     for (i = 0 ; i < n ; i++)
88     {
89 	Rp [i+1] = Rp [i] + W [i] ;
90     }
91     for (i = 0 ; i < n ; i++)
92     {
93 	W [i] = Rp [i] ;
94 	Flag [i] = EMPTY ;
95     }
96 
97     /* --------------------------------------------------------------------- */
98     /* construct the row form matrix R */
99     /* --------------------------------------------------------------------- */
100 
101     /* R = row form of pattern of A */
102     for (j = 0 ; j < n ; j++)
103     {
104 	p2 = Ap [j+1] ;
105 	for (p = Ap [j] ; p < p2 ; p++)
106 	{
107 	    i = Ai [p] ;
108 	    if (Flag [i] != j)
109 	    {
110 		/* row index i has not yet appeared in column j */
111 		Ri [W [i]++] = j ;  /* put col j in row i */
112 		Flag [i] = j ;	    /* flag row index i as appearing in col j*/
113 	    }
114 	}
115     }
116 
117 version(NDEBUG){}
118 else {
119     ASSERT (AMD_valid (n, n, Rp, Ri) == AMD_OK) ;
120     for (j = 0 ; j < n ; j++)
121     {
122 	ASSERT (W [j] == Rp [j+1]) ;
123     }
124 }
125 }