cs_maxtrans.c
Go to the documentation of this file.
00001 #include "cs.h"
00002 /* find an augmenting path starting at column k and extend the match if found */
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 ;                        /* start with just node k in jstack */
00008     while (head >= 0)
00009     {
00010         /* --- Start (or continue) depth-first-search at node j ------------- */
00011         j = js [head] ;                 /* get j from top of jstack */
00012         if (w [j] != k)                 /* 1st time j visited for kth path */
00013         {
00014             w [j] = k ;                 /* mark j as visited for kth path */
00015             for (p = cheap [j] ; p < Ap [j+1] && !found ; p++)
00016             {
00017                 i = Ai [p] ;            /* try a cheap assignment (i,j) */
00018                 found = (jmatch [i] == -1) ;
00019             }
00020             cheap [j] = p ;             /* start here next time j is traversed*/
00021             if (found)
00022             {
00023                 is [head] = i ;         /* column j matched with row i */
00024                 break ;                 /* end of augmenting path */
00025             }
00026             ps [head] = Ap [j] ;        /* no cheap match: start dfs for j */
00027         }
00028         /* --- Depth-first-search of neighbors of j ------------------------- */
00029         for (p = ps [head] ; p < Ap [j+1] ; p++)
00030         {
00031             i = Ai [p] ;                /* consider row i */
00032             if (w [jmatch [i]] == k) continue ; /* skip jmatch [i] if marked */
00033             ps [head] = p + 1 ;         /* pause dfs of node j */
00034             is [head] = i ;             /* i will be matched with j if found */
00035             js [++head] = jmatch [i] ;  /* start dfs at column jmatch [i] */
00036             break ;
00037         }
00038         if (p == Ap [j+1]) head-- ;     /* node j is done; pop from stack */
00039     }                                   /* augment the match if path found: */
00040     if (found) for (p = head ; p >= 0 ; p--) jmatch [is [p]] = js [p] ;
00041 }
00042 
00043 /* find a maximum transveral */
00044 int *cs_maxtrans (const cs *A, int seed)  /*[jmatch [0..m-1]; imatch [0..n-1]]*/
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) ;                /* check inputs */
00050     n = A->n ; m = A->m ; Ap = A->p ; Ai = A->i ;
00051     w = jimatch = (int*)cs_calloc (m+n, sizeof (int)) ;   /* allocate result */
00052     if (!jimatch) return (NULL) ;
00053     for (k = 0, j = 0 ; j < n ; j++)    /* count nonempty rows and columns */
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]) ;        /* count entries already on diagonal */
00060         }
00061     }
00062     if (k == CS_MIN (m,n))              /* quick return if diagonal zero-free */
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) ; /* transpose if needed */
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 = (int*)cs_malloc (5*n, sizeof (int)) ;             /* get workspace */
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] ;  /* for cheap assignment */
00081     for (j = 0 ; j < n ; j++) w [j] = -1 ;          /* all columns unflagged */
00082     for (i = 0 ; i < m ; i++) jmatch [i] = -1 ;     /* nothing matched yet */
00083     q = cs_randperm (n, seed) ;                     /* q = random permutation */
00084     for (k = 0 ; k < n ; k++)   /* augment, starting at column q[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 ;     /* find row match */
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 }


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:36:56