00001 #include "cs.h"
00002
00003 static void cs_augment (int k, const cs *A, int *jmatch, int *cheap, int *w,
00004 int *js, int *is, int *ps)
00005 {
00006 int found = 0, p, i = -1, *Ap = A->p, *Ai = A->i, head = 0, j ;
00007 js [0] = k ;
00008 while (head >= 0)
00009 {
00010
00011 j = js [head] ;
00012 if (w [j] != k)
00013 {
00014 w [j] = k ;
00015 for (p = cheap [j] ; p < Ap [j+1] && !found ; p++)
00016 {
00017 i = Ai [p] ;
00018 found = (jmatch [i] == -1) ;
00019 }
00020 cheap [j] = p ;
00021 if (found)
00022 {
00023 is [head] = i ;
00024 break ;
00025 }
00026 ps [head] = Ap [j] ;
00027 }
00028
00029 for (p = ps [head] ; p < Ap [j+1] ; p++)
00030 {
00031 i = Ai [p] ;
00032 if (w [jmatch [i]] == k) continue ;
00033 ps [head] = p + 1 ;
00034 is [head] = i ;
00035 js [++head] = jmatch [i] ;
00036 break ;
00037 }
00038 if (p == Ap [j+1]) head-- ;
00039 }
00040 if (found) for (p = head ; p >= 0 ; p--) jmatch [is [p]] = js [p] ;
00041 }
00042
00043
00044 int *cs_maxtrans (const cs *A, int seed)
00045 {
00046 int i, j, k, n, m, p, n2 = 0, m2 = 0, *Ap, *jimatch, *w, *cheap, *js, *is,
00047 *ps, *Ai, *Cp, *jmatch, *imatch, *q ;
00048 cs *C ;
00049 if (!CS_CSC (A)) return (NULL) ;
00050 n = A->n ; m = A->m ; Ap = A->p ; Ai = A->i ;
00051 w = jimatch = cs_calloc (m+n, sizeof (int)) ;
00052 if (!jimatch) return (NULL) ;
00053 for (k = 0, j = 0 ; j < n ; j++)
00054 {
00055 n2 += (Ap [j] < Ap [j+1]) ;
00056 for (p = Ap [j] ; p < Ap [j+1] ; p++)
00057 {
00058 w [Ai [p]] = 1 ;
00059 k += (j == Ai [p]) ;
00060 }
00061 }
00062 if (k == CS_MIN (m,n))
00063 {
00064 jmatch = jimatch ; imatch = jimatch + m ;
00065 for (i = 0 ; i < k ; i++) jmatch [i] = i ;
00066 for ( ; i < m ; i++) jmatch [i] = -1 ;
00067 for (j = 0 ; j < k ; j++) imatch [j] = j ;
00068 for ( ; j < n ; j++) imatch [j] = -1 ;
00069 return (cs_idone (jimatch, NULL, NULL, 1)) ;
00070 }
00071 for (i = 0 ; i < m ; i++) m2 += w [i] ;
00072 C = (m2 < n2) ? cs_transpose (A,0) : ((cs *) A) ;
00073 if (!C) return (cs_idone (jimatch, (m2 < n2) ? C : NULL, NULL, 0)) ;
00074 n = C->n ; m = C->m ; Cp = C->p ;
00075 jmatch = (m2 < n2) ? jimatch + n : jimatch ;
00076 imatch = (m2 < n2) ? jimatch : jimatch + m ;
00077 w = cs_malloc (5*n, sizeof (int)) ;
00078 if (!w) return (cs_idone (jimatch, (m2 < n2) ? C : NULL, w, 0)) ;
00079 cheap = w + n ; js = w + 2*n ; is = w + 3*n ; ps = w + 4*n ;
00080 for (j = 0 ; j < n ; j++) cheap [j] = Cp [j] ;
00081 for (j = 0 ; j < n ; j++) w [j] = -1 ;
00082 for (i = 0 ; i < m ; i++) jmatch [i] = -1 ;
00083 q = cs_randperm (n, seed) ;
00084 for (k = 0 ; k < n ; k++)
00085 {
00086 cs_augment (q ? q [k]: k, C, jmatch, cheap, w, js, is, ps) ;
00087 }
00088 cs_free (q) ;
00089 for (j = 0 ; j < n ; j++) imatch [j] = -1 ;
00090 for (i = 0 ; i < m ; i++) if (jmatch [i] >= 0) imatch [jmatch [i]] = i ;
00091 return (cs_idone (jimatch, (m2 < n2) ? C : NULL, w, 1)) ;
00092 }