← Back to branch summary

~logan/ubuntu/trusty/suitesparse/4.2.1-3ubuntu1

« back to all changes in this revision

Viewing changes to CXSparse/MATLAB/CSparse/cs_lsolve_mex.c

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2007-05-29 09:36:29 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070529093629-zowquo0b7slkk6nc
Tags: 3.0.0-2
http://bugs.debian.org/426349
http://bugs.debian.org/426354
* suitesparse builds properly twice in a row
* Bug fix: "suitesparse - FTBFS: Broken build depens: libgfortran1-dev",
  thanks to Bastian Blank (Closes: #426349).
* Bug fix: "suitesparse_3.0.0-1: FTBFS: build-depends on
  libgfortran1-dev", thanks to Steve Langasek (Closes: #426354).

Show diffs side-by-side

added added

removed removed

Lines of Context:
CXSparse/MATLAB/CSparse/cs_lsolve_mex.c
 
1
#include "cs_mex.h"
 
2
/* cs_lsolve: x=L\b.  L must be sparse and lower triangular.  b must be a
 
3
 * full or sparse vector.  x is full or sparse, depending on b.
 
4
 *
 
5
 * Time taken is O(flop count), which may be less than n if b is sparse,
 
6
 * depending on L and b.
 
7
 *
 
8
 * This function works with MATLAB 7.2, but is not perfectly compatible with
 
9
 * the requirements of a MATLAB mexFunction when b is sparse.  X is returned
 
10
 * as an unsorted sparse vector.  Also, this mexFunction temporarily modifies
 
11
 * its input, L, by modifying L->p (in the cs_dfs function) and then restoring
 
12
 * it.  This could be corrected by creating a copy of L->p
 
13
 * (see cs_dmperm_mex.c), but this would take O(n) time, destroying the
 
14
 * O(flop count) time complexity of this function.
 
15
 *
 
16
 * Note that b cannot be sparse complex.  This function does not support
 
17
 * sparse complex L and b because the sparse x=L\b only accesses part of the
 
18
 * matrix L.  Converting L from a MATLAB complex matrix to a CXSparse complex
 
19
 * matrix requires all of L to be accessed, defeating the purpose of this
 
20
 * function.
 
21
 *
 
22
 * L can be sparse complex, but in that case b must be full real or complex,
 
23
 * not sparse.
 
24
 */
 
25
 
 
26
void mexFunction
 
27
(
 
28
    int nargout,
 
29
    mxArray *pargout [ ],
 
30
    int nargin,
 
31
    const mxArray *pargin [ ]
 
32
)
 
33
{
 
34
    CS_INT top, nz, p, *xi, n ;
 
35
    if (nargout > 1 || nargin != 2)
 
36
    {
 
37
        mexErrMsgTxt ("Usage: x = cs_lsolve(L,b)") ;
 
38
    }
 
39
    if (mxIsSparse (pargin [1]))
 
40
    {
 
41
        cs_dl Lmatrix, Bmatrix, *L, *B, *X ;
 
42
        double *x ;
 
43
        if (mxIsComplex (pargin [0]) || mxIsComplex (pargin [1]))
 
44
        {
 
45
            mexErrMsgTxt ("sparse complex case not supported") ;
 
46
        }
 
47
        L = cs_dl_mex_get_sparse (&Lmatrix, 1, 1, pargin [0]) ;/* get L */
 
48
        n = L->n ;
 
49
        B = cs_dl_mex_get_sparse (&Bmatrix, 0, 1, pargin [1]) ;/* get sparse b*/
 
50
        cs_mex_check (0, n, 1, 0, 1, 1, pargin [1]) ;
 
51
        xi = cs_dl_malloc (2*n, sizeof (CS_INT)) ;          /* get workspace */
 
52
        x  = cs_dl_malloc (n, sizeof (double)) ;
 
53
        top = cs_dl_spsolve (L, B, 0, xi, x, NULL, 1) ;     /* x = L\b */
 
54
        X = cs_dl_spalloc (n, 1, n-top, 1, 0) ;     /* create sparse x*/
 
55
        X->p [0] = 0 ;
 
56
        nz = 0 ;
 
57
        for (p = top ; p < n ; p++)
 
58
        {
 
59
            X->i [nz] = xi [p] ;
 
60
            X->x [nz++] = x [xi [p]] ;
 
61
        }
 
62
        X->p [1] = nz ;
 
63
        pargout [0] = cs_dl_mex_put_sparse (&X) ;
 
64
        cs_free (x) ;
 
65
        cs_free (xi) ;
 
66
    }
 
67
    else if (mxIsComplex (pargin [0]) || mxIsComplex (pargin [1]))
 
68
    {
 
69
#ifndef NCOMPLEX
 
70
        cs_cl Lmatrix, Bmatrix, *L, *B, *X ;
 
71
        cs_complex_t *x ;
 
72
        L = cs_cl_mex_get_sparse (&Lmatrix, 1, pargin [0]) ;    /* get L */
 
73
        n = L->n ;
 
74
        x = cs_cl_mex_get_double (n, pargin [1]) ;              /* x = b */
 
75
        cs_cl_lsolve (L, x) ;                                   /* x = L\x */
 
76
        cs_free (L->x) ;
 
77
        pargout [0] = cs_cl_mex_put_double (n, x) ;             /* return x */
 
78
#else
 
79
        mexErrMsgTxt ("complex matrices not supported") ;
 
80
#endif
 
81
    }
 
82
    else
 
83
    {
 
84
        cs_dl Lmatrix, *L ;
 
85
        double *x, *b ;
 
86
        L = cs_dl_mex_get_sparse (&Lmatrix, 1, 1, pargin [0]) ; /* get L */
 
87
        n = L->n ;
 
88
        b = cs_dl_mex_get_double (n, pargin [1]) ;              /* get b */
 
89
        x = cs_dl_mex_put_double (n, b, &(pargout [0])) ;       /* x = b */
 
90
        cs_dl_lsolve (L, x) ;                                   /* x = L\x */
 
91
    }
 
92
}