← Back to branch summary

~ubuntu-branches/ubuntu/natty/suitesparse/natty

« back to all changes in this revision

Viewing changes to KLU/MATLAB/old/kluzmex.c

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2006-12-22 10:16:15 UTC
  • Revision ID: james.westby@ubuntu.com-20061222101615-2ohaj8902oix2rnk
Tags: upstream-2.3.1
Import upstream version 2.3.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
KLU/MATLAB/old/kluzmex.c
 
1
/* ========================================================================== */
 
2
/* === klu mexFunction ====================================================== */
 
3
/* ========================================================================== */
 
4
 
 
5
/* Factorizes A(p,q) into L*U.  Usage:
 
6
 *
 
7
 * [L,U,p,t] = klu (A)                  % Q assumed to be identity
 
8
 * [L,U,p,t] = klu (A, q)
 
9
 * [L,U,p,t] = klu (A, q, Control)
 
10
 *
 
11
 * Control (1): threshold partial pivoting tolerance (0 means force diagonal
 
12
 *              pivoting, 1 means conventional partial pivoting.  A value of
 
13
 *              0.5 means the diagonal will be picked if its absolute value
 
14
 *              is >= 0.5 times the largest value in the column; otherwise
 
15
 *              the largest value is selected.  Default: 1.0.
 
16
 * Control (2): if positive, this is the initial size of the L matrix, in # of
 
17
 *              nonzero entries.  If negative, the initial size of L is
 
18
 *              (-Control (2) * nnz (A)).  Default: -10.
 
19
 * Control (3): if positive, this is the initial size of the U matrix, in # of
 
20
 *              nonzero entries.  If negative, the initial size of U is
 
21
 *              (-Control (2) * nnz (A)).  Default: -10.
 
22
 * Control (4): memory growth factor
 
23
 *
 
24
 * t = [cputime noffdiag umin umax], output statistics.
 
25
 *
 
26
 * If Control is not present, or present but not of the correct size,
 
27
 * defaults are used.
 
28
 *
 
29
 * If klu is compiled with the default options (NRECIPROCAL not defined), then
 
30
 * the diagonal of U is returned inverted on output.
 
31
 */
 
32
 
 
33
/* ========================================================================== */
 
34
 
 
35
/* #include "klu.h" */
 
36
#include "klu_internal.h"
 
37
#include "mex.h"
 
38
#include "tictoc.h"
 
39
 
 
40
void mexFunction
 
41
(
 
42
    int nargout,
 
43
    mxArray *pargout [ ],
 
44
    int nargin,
 
45
    const mxArray *pargin [ ]
 
46
)
 
47
{
 
48
    double Control [KLU_CONTROL], tt [2], umin, umax ;
 
49
    double *Ax, *Px, *User_Control, *Lx, *Ux, *T, *Qx ;
 
50
    _Complex double *Atemp, *LU, *X ;
 
51
    int n, result, k, anz, s, lnz, unz, j, p, col, noffdiag, nrealloc ;
 
52
    int *Ap, *Ai, *Li, *Ui, *P, *Pp, *Pi, 
 
53
        *Q, *Work, *Llen, *Ulen, *Lip, *Uip, sing_col ;
 
54
 
 
55
    /* ---------------------------------------------------------------------- */
 
56
    /* get inputs */
 
57
    /* ---------------------------------------------------------------------- */
 
58
 
 
59
    if (nargin < 1 || nargin > 3 || !(nargout == 3 || nargout == 4))
 
60
    {
 
61
        mexErrMsgTxt ("Usage: [L,U,P,t] = klu (A,Q,Control), where t, Q, and Control are optional.") ;
 
62
    }
 
63
    n = mxGetM (pargin [0]) ;
 
64
    if (!mxIsSparse (pargin [0]) || n != mxGetN (pargin [0])
 
65
        || n == 0)
 
66
    {
 
67
        mexErrMsgTxt ("klu: A must be sparse, square, complex, and non-empty") ;
 
68
    }
 
69
 
 
70
    /* get sparse matrix A */
 
71
    Ap = mxGetJc (pargin [0]) ;
 
72
    Ai = mxGetIr (pargin [0]) ;
 
73
    Ax = mxGetPr (pargin [0]) ;
 
74
    Az = mxGetPi (pargin [0]) ;
 
75
    anz = Ap [n] ;
 
76
    Atemp = mxMalloc(anz * sizeof(_Complex double)) ;
 
77
    if (!Atemp)             
 
78
        mexErrMsgTxt("malloc failed") ;
 
79
    for (k = 0; k < anz ; k++)
 
80
        Atemp[k] = Ax [k] + Az [k] * _Complex_I ;
 
81
    /* get input column permutation Q */
 
82
    if (nargin > 1)
 
83
    {
 
84
        if (!mxIsDouble (pargin [1]) || n != mxGetNumberOfElements (pargin [1]))
 
85
        {
 
86
            mexErrMsgTxt ("klu: Q must be a dense 1-by-n vector") ; 
 
87
        }
 
88
        Qx = mxGetPr (pargin [1]) ;
 
89
        Q = (int *) mxMalloc (n * sizeof (int)) ;
 
90
        for (k = 0 ; k < n ; k++)
 
91
        {
 
92
            col = (int) (Qx [k]) - 1 ;  /* convert from 1-based to 0-based */
 
93
            if (col < 0 || col >= n)
 
94
            {
 
95
                mexErrMsgTxt ("klu: Q not a valid permutation\n") ;
 
96
            }
 
97
            Q [k] = col ;
 
98
        }
 
99
    }
 
100
    else
 
101
    {
 
102
        /* klu will assume that Q is the identity permutation */
 
103
        Q = (int *) NULL ;
 
104
    }
 
105
 
 
106
    /* get control parameters */
 
107
    klu_defaults (Control) ;
 
108
    if (nargin > 2)
 
109
    {
 
110
        if (!mxIsDouble (pargin [2]))
 
111
        {
 
112
            mexErrMsgTxt ("klu: Control must be real") ;
 
113
        }
 
114
        User_Control = mxGetPr (pargin [2]) ;
 
115
        s = mxGetNumberOfElements (pargin [2]) ;
 
116
        for (k = 0 ; k < s ; k++)
 
117
        {
 
118
            Control [k] = User_Control [k] ;
 
119
        }
 
120
    }
 
121
 
 
122
    P  = (int *) mxMalloc (n * sizeof (int)) ;
 
123
    Llen = (int *) mxMalloc (n * sizeof (int)) ;
 
124
    Ulen = (int *) mxMalloc (n * sizeof (int)) ;
 
125
    Lip = (int *) mxMalloc ((n+1) * sizeof (int)) ;
 
126
    Uip = (int *) mxMalloc ((n+1) * sizeof (int)) ;
 
127
 
 
128
    X    = (_Complex double *) mxMalloc (n * sizeof (_Complex double)) ;
 
129
    Work = (int *) mxMalloc (5*n * sizeof (int)) ;
 
130
 
 
131
    /* ---------------------------------------------------------------------- */
 
132
    /* factorize */
 
133
    /* ---------------------------------------------------------------------- */
 
134
 
 
135
    /* my_tic (tt) ; */
 
136
 
 
137
    result = klu_factor (n, Ap, Ai, Atemp, Q, (double *) NULL,
 
138
                 &LU, Llen, Ulen, Lip, Uip, P, &lnz, &unz, 
 
139
                 &noffdiag, &umin, &umax, &nrealloc, &sing_col, X, Work, 
 
140
                 /* no BTF or scaling here */
 
141
                 0, (int *) NULL, (double *) NULL, 0, 0, (int *) NULL, (int *) NULL,
 
142
                 (double *) NULL) ;
 
143
 
 
144
    /* my_toc (tt) ; */
 
145
 
 
146
    mxFree (X) ;
 
147
    mxFree (Work) ;
 
148
 
 
149
    if (nargout == 4)
 
150
    {
 
151
        pargout [3] = mxCreateDoubleMatrix (1, 4, mxREAL) ;
 
152
        T = mxGetPr (pargout [3]) ;
 
153
        T [0] = tt [1] ;
 
154
        T [1] = noffdiag ;
 
155
        T [2] = umin ;
 
156
        T [3] = umax ;
 
157
    }
 
158
 
 
159
    if (result == KLU_OUT_OF_MEMORY)
 
160
    {
 
161
        mexErrMsgTxt ("klu: out of memory") ;
 
162
    }
 
163
 
 
164
    /* NOTE: this should be done in a better way, without copying, but when
 
165
     * I tried to do it I got segmentation faults, and gave up ... */
 
166
 
 
167
    /* create sparse matrix for L */
 
168
    pargout [0] = mxCreateSparse (n, n, lnz, mxREAL) ;
 
169
    mxSetJc (pargout [0], Lip);
 
170
    Li = mxGetIr (pargout [0]) ;
 
171
    Lx = mxGetPr (pargout [0]) ;
 
172
    KLU_z_convert (LU, Lip, Llen, Li, Lx, n) ; 
 
173
 
 
174
    /* create sparse matrix for U */
 
175
    pargout [1] = mxCreateSparse (n, n, unz, mxREAL) ;
 
176
    mxSetJc (pargout [1], Uip);
 
177
    Ui = mxGetIr (pargout [1]) ;
 
178
    Ux = mxGetPr (pargout [1]) ;
 
179
    KLU_z_convert (LU, Uip, Ulen, Ui, Ux, n) ; 
 
180
 
 
181
    mxFree (LU) ; 
 
182
    mxFree (Llen) ; 
 
183
    mxFree (Ulen) ; 
 
184
 
 
185
    /* create permutation vector for P */
 
186
    pargout [2] = mxCreateDoubleMatrix (1, n, mxREAL) ;
 
187
    Px = mxGetPr (pargout [2]) ;
 
188
    for (k = 0 ; k < n ; k++)
 
189
    {
 
190
        Px [k] = P [k] + 1 ;    /* convert to 1-based */
 
191
    }
 
192
 
 
193
    mxFree (P) ;
 
194
 
 
195
}