« back to all changes in this revision
Viewing changes to SPQR/Source/SuiteSparseQR_expert.cpp
- browse files at revision 8
- compare with another revision
- download diff
- download tarball
- view history from revision 8
- Committer: Bazaar Package Importer
- Author(s): Rafael Laboissiere
- Date: 2009-02-24 11:08:12 UTC
- mfrom: (7.2.2 sid)
- Revision ID: james.westby@ubuntu.com-20090224110812-hawvr3dh5kjbvlae
debian/control: Add an epoch to the version number of
libsuitesparse-3.0.2 in replaces/conflicts for libcolamd-3.2.0
(really, closes: #516725)
libsuitesparse-3.0.2 in replaces/conflicts for libcolamd-3.2.0
(really, closes: #516725)
- files added:
- MATLAB_Tools/GEE
- MATLAB_Tools/GEE/private
- MATLAB_Tools/find_components
- MATLAB_Tools/spok
- MATLAB_Tools/spok/private
- SPQR
- SPQR/Demo
- SPQR/Doc
- SPQR/Include
- SPQR/Lib
- SPQR/MATLAB
- SPQR/Matrix
- SPQR/Source
- SPQR/Tcov
- files removed:
- debian/libsuitesparse
- debian/libsuitesparse/usr
- debian/libsuitesparse/usr/share
- debian/libsuitesparse/usr/share/lintian
- debian/libsuitesparse/usr/share/lintian/overrides
- files modified:
- CHOLMOD/Check/cholmod_check.c
added
removed
SPQR/Source/SuiteSparseQR_expert.cpp
1
// =============================================================================
2
// === SuiteSparseQR_expert ====================================================
3
// =============================================================================
4
5
// This file includes a suite of "expert" user-callable routines that allow the
6
// reuse of a QR factorization, including the ability to find the min 2-norm
7
// solution of an under-determined system of equations:
8
//
9
// The QR factorization can be computed in one of two ways:
10
//
11
// SuiteSparseQR_factorize QR factorization, returning a QR object
12
// This performs both symbolic analysis and
13
// numeric factorization, and exploits singletons.
14
//
15
// SuiteSparseQR_symbolic symbolic QR factorization, based purely on the
16
// nonzero pattern of A; to be followed by:
17
// SuiteSparseQR_numeric numeric QR factorization. Does not exploit
18
// singletons. Note that H is always kept.
19
//
20
// SuiteSparseQR_solve forward/backsolve using R from the QR object
21
// SuiteSparseQR_qmult multiply by Q or Q', using Q from the QR object
22
// SuiteSparseQR_min2norm min 2-norm solution for x=A\b
23
// SuiteSparseQR_free free the QR object
24
//
25
// If all you need is to solve a least-squares problem, or to compute a QR
26
// factorization and return the results in vanilla sparse matrix format, see
27
// the routines in SuiteSparseQR.cpp and SuiteSparseQR_qmult.cpp instead.
28
//
29
// NOTE: Q the matrix is not formed; H the Householder vectors are always kept.
30
31
#ifndef NEXPERT
32
#include "spqr.hpp"
33
34
// =============================================================================
35
// === SuiteSparseQR_symbolic ==================================================
36
// =============================================================================
37
38
// This returns a QR factorization object with a NULL numeric part. It must
39
// be followed by a numeric factorization, by SuiteSparseQR_numeric.
40
41
template <typename Entry>
42
SuiteSparseQR_factorization <Entry> *SuiteSparseQR_symbolic
43
(
44
// inputs:
45
int ordering, // all, except 3:given treated as 0:fixed
46
int allow_tol, // if FALSE, tol is ignored by the numeric
47
// factorization, and no rank detection is performed
48
cholmod_sparse *A, // sparse matrix to factorize (A->x ignored)
49
cholmod_common *cc // workspace and parameters
50
)
51
{
52
#ifdef TIMING
53
double t0 = spqr_time ( ) ;
54
#endif
55
RETURN_IF_NULL_COMMON (NULL) ;
56
RETURN_IF_NULL (A, NULL) ;
57
Int xtype = spqr_type <Entry> ( ) ;
58
RETURN_IF_XTYPE_INVALID (A, NULL) ;
59
cc->status = CHOLMOD_OK ;
60
61
// -------------------------------------------------------------------------
62
// get inputs and allocate result
63
// -------------------------------------------------------------------------
64
65
SuiteSparseQR_factorization <Entry> *QR ;
66
spqr_symbolic *QRsym ;
67
68
QR = (SuiteSparseQR_factorization <Entry> *)
69
cholmod_l_malloc (1, sizeof (SuiteSparseQR_factorization <Entry>), cc) ;
70
71
if (cc->status < CHOLMOD_OK)
72
{
73
// out of memory
74
return (NULL) ;
75
}
76
77
// -------------------------------------------------------------------------
78
// perform the symbolic analysis
79
// -------------------------------------------------------------------------
80
81
QR->QRsym = QRsym = spqr_analyze (A, ordering, NULL, allow_tol, TRUE, cc) ;
82
83
QR->QRnum = NULL ; // allocated later, by numeric factorization
84
85
// singleton information:
86
QR->R1p = NULL ; // no singletons; these always remain NULL or 0
87
QR->R1j = NULL ;
88
QR->R1x = NULL ;
89
QR->P1inv = NULL ;
90
QR->HP1inv = NULL ;
91
QR->r1nz = 0 ;
92
QR->n1rows = 0 ;
93
QR->n1cols = 0 ;
94
cc->SPQR_istat [5] = 0 ; // number of columns singletons
95
cc->SPQR_istat [6] = 0 ; // number of singleton rows
96
97
QR->Q1fill = NULL ; // may change below
98
99
QR->Rmap = NULL ; // may be allocated by numeric factorization
100
QR->RmapInv = NULL ;
101
102
QR->narows = A->nrow ;
103
QR->nacols = A->ncol ;
104
QR->bncols = 0 ; // [A B] is not factorized
105
106
QR->allow_tol = (allow_tol != 0) ;
107
QR->tol = QR->allow_tol ? SPQR_DEFAULT_TOL : EMPTY ;
108
109
if (cc->status < CHOLMOD_OK)
110
{
111
// out of memory
112
spqr_freefac (&QR, cc) ;
113
return (NULL) ;
114
}
115
116
// -------------------------------------------------------------------------
117
// copy the fill-reducing ordering
118
// -------------------------------------------------------------------------
119
120
if (QRsym->Qfill != NULL)
121
{
122
Int n, k, *Qfill, *Q1fill ;
123
Qfill = QRsym->Qfill ;
124
n = A->ncol ;
125
Q1fill = (Int *) cholmod_l_malloc (n, sizeof (Int), cc) ;
126
QR->Q1fill = Q1fill ;
127
if (cc->status < CHOLMOD_OK)
128
{
129
// out of memory
130
spqr_freefac (&QR, cc) ;
131
return (NULL) ;
132
}
133
for (k = 0 ; k < n ; k++)
134
{
135
Q1fill [k] = Qfill [k] ;
136
}
137
}
138
139
#ifdef TIMING
140
double t1 = spqr_time ( ) ;
141
cc->other1 [1] = t1 - t0 ;
142
#endif
143
144
return (QR) ;
145
}
146
147
template
148
SuiteSparseQR_factorization <double> *SuiteSparseQR_symbolic <double>
149
(
150
// inputs:
151
int ordering, // all, except 3:given treated as 0:fixed
152
int allow_tol, // if FALSE, tol is ignored by the numeric
153
// factorization, and no rank detection is performed
154
cholmod_sparse *A, // sparse matrix to factorize (A->x ignored)
155
cholmod_common *cc // workspace and parameters
156
) ;
157
158
template
159
SuiteSparseQR_factorization <Complex> *SuiteSparseQR_symbolic <Complex>
160
(
161
// inputs:
162
int ordering, // all, except 3:given treated as 0:fixed
163
int allow_tol, // if FALSE, tol is ignored by the numeric
164
// factorization, and no rank detection is performed
165
cholmod_sparse *A, // sparse matrix to factorize (A->x ignored)
166
cholmod_common *cc // workspace and parameters
167
) ;
168
169
// =============================================================================
170
// === SuiteSparseQR_numeric ===================================================
171
// =============================================================================
172
173
// numeric QR factorization; must be preceded by a call to
174
// SuiteSparseQR_symbolic. Alternatively, if SuiteSparseQR_factorize
175
// did not find any singletons (QR->n1cols == 0), a call to
176
// SuiteSparseQR_factorize can be followed by a factorization by this function.
177
178
template <typename Entry> int SuiteSparseQR_numeric
179
(
180
// inputs:
181
double tol, // treat columns with 2-norm <= tol as zero
182
cholmod_sparse *A, // sparse matrix to factorize
183
// input/output
184
SuiteSparseQR_factorization <Entry> *QR,
185
cholmod_common *cc // workspace and parameters
186
)
187
{
188
#ifdef TIMING
189
double t0 = spqr_time ( ) ;
190
#endif
191
RETURN_IF_NULL_COMMON (FALSE) ;
192
RETURN_IF_NULL (A, FALSE) ;
193
RETURN_IF_NULL (QR, FALSE) ;
194
Int xtype = spqr_type <Entry> ( ) ;
195
RETURN_IF_XTYPE_INVALID (A, FALSE) ;
196
cc->status = CHOLMOD_OK ;
197
198
if (QR->n1cols > 0 || QR->bncols > 0)
199
{
200
// this numeric factorization phase cannot be used if singletons were
201
// found, or if [A B] was factorized.
202
ERROR (CHOLMOD_INVALID, "cannot refactorize w/singletons or [A B]") ;
203
return (FALSE) ;
204
}
205
206
Int n = A->ncol ;
207
208
// -------------------------------------------------------------------------
209
// get the column 2-norm tolerance
210
// -------------------------------------------------------------------------
211
212
if (QR->allow_tol)
213
{
214
// compute default tol, if requested
215
if (tol <= SPQR_DEFAULT_TOL)
216
{
217
tol = spqr_tol <Entry> (A, cc) ;
218
}
219
}
220
else
221
{
222
// tol is ignored; no rank detection performed
223
tol = EMPTY ;
224
}
225
QR->tol = tol ;
226
227
// -------------------------------------------------------------------------
228
// numeric factorization
229
// -------------------------------------------------------------------------
230
231
// free the existing numeric factorization, if any
232
spqr_freenum (&(QR->QRnum), cc) ;
233
234
// compute the new factorization
235
QR->QRnum = spqr_factorize <Entry> (&A, FALSE, tol, n, QR->QRsym, cc) ;
236
237
if (cc->status < CHOLMOD_OK)
238
{
239
// out of memory; QR factorization remains a symbolic-only object
240
ASSERT (QR->QRnum == NULL) ;
241
return (FALSE) ;
242
}
243
244
QR->rank = QR->QRnum->rank1 ;
245
246
// -------------------------------------------------------------------------
247
// find the mapping for the squeezed R, if A is rank deficient
248
// -------------------------------------------------------------------------
249
250
if (QR->rank < n && !spqr_rmap (QR, cc))
251
{
252
// out of memory; QR factorization remains a symbolic-only object
253
spqr_freenum (&(QR->QRnum), cc) ;
254
return (FALSE) ;
255
}
256
257
// -------------------------------------------------------------------------
258
// output statistics
259
// -------------------------------------------------------------------------
260
261
cc->SPQR_istat [4] = QR->rank ; // estimated rank of A
262
cc->SPQR_xstat [1] = tol ; // tol used
263
264
#ifdef TIMING
265
double t1 = spqr_time ( ) ;
266
cc->other1 [2] = t1 - t0 ;
267
#endif
268
269
return (TRUE) ;
270
}
271
272
template int SuiteSparseQR_numeric <double>
273
(
274
// inputs:
275
double tol, // treat columns with 2-norm <= tol as zero
276
cholmod_sparse *A, // sparse matrix to factorize
277
// input/output
278
SuiteSparseQR_factorization <double> *QR,
279
cholmod_common *cc // workspace and parameters
280
) ;
281
282
template int SuiteSparseQR_numeric <Complex>
283
(
284
// inputs:
285
double tol, // treat columns with 2-norm <= tol as zero
286
cholmod_sparse *A, // sparse matrix to factorize
287
// input/output
288
SuiteSparseQR_factorization <Complex> *QR,
289
cholmod_common *cc // workspace and parameters
290
) ;
291
292
// =============================================================================
293
// === SuiteSparseQR_factorize =================================================
294
// =============================================================================
295
296
// QR factorization of a sparse matrix A, including symbolic ordering/analysis,
297
// and numeric factorization. This function is a combined symbolic+numeric
298
// factorization, because it exploits singletons. It factorizes A*E into Q*R,
299
// where E is the fill-reducing ordering, Q is in Householder form, and R is
300
// composed of a set of singleton rows followed by the rows of the multifrontal
301
// factorization of the singleton-free submatrix.
302
//
303
// The MATLAB equivalent of this function is [Q,R,E]=spqr(A,opts) where opts.Q
304
// = 'Householder', except that Q, R, and E are all stored in a single QR
305
// object (of type SuiteSparseQR_factorization), rather than being put into
306
// MATLAB format. That extraction is done by the overloaded SuiteSparseQR
307
// function, which also frees its QR object after extracting Q, R, and E.
308
//
309
// If you have a least-squares problem with a single right-hand-side B (where B
310
// is a vector or matrix, or sparse or dense), and do not need to reuse the QR
311
// factorization of A, then use the SuiteSparseQR function instead (which
312
// corresponds to the MATLAB statement x=A\B, using a sparse QR factorization).
313
314
template <typename Entry>
315
SuiteSparseQR_factorization <Entry> *SuiteSparseQR_factorize
316
(
317
// inputs, not modified:
318
int ordering, // all, except 3:given treated as 0:fixed
319
double tol, // treat columns with 2-norm <= tol as zero
320
cholmod_sparse *A, // sparse matrix to factorize
321
cholmod_common *cc // workspace and parameters
322
)
323
{
324
RETURN_IF_NULL_COMMON (NULL) ;
325
RETURN_IF_NULL (A, NULL) ;
326
Int xtype = spqr_type <Entry> ( ) ;
327
RETURN_IF_XTYPE_INVALID (A, NULL) ;
328
cc->status = CHOLMOD_OK ;
329
// B is not present, and always keep H:
330
return (spqr_1factor <Entry> (ordering, tol, 0, TRUE, A,
331
0, NULL, NULL, NULL, cc)) ;
332
}
333
334
template SuiteSparseQR_factorization <double> *SuiteSparseQR_factorize <double>
335
(
336
// inputs, not modified:
337
int ordering, // all, except 3:given treated as 0:fixed
338
double tol, // treat columns with 2-norm <= tol as zero
339
cholmod_sparse *A, // sparse matrix to factorize
340
// workspace and parameters
341
cholmod_common *cc
342
) ;
343
344
template SuiteSparseQR_factorization <Complex> *SuiteSparseQR_factorize<Complex>
345
(
346
// inputs, not modified:
347
int ordering, // all, except 3:given treated as 0:fixed
348
double tol, // treat columns with 2-norm <= tol as zero
349
cholmod_sparse *A, // sparse matrix to factorize
350
// workspace and parameters
351
cholmod_common *cc
352
) ;
353
354
355
// =============================================================================
356
// === rtsolve =================================================================
357
// =============================================================================
358
359
// Solve X = R'\(E'*B) or X = R'\B using the QR factorization from spqr_1factor
360
// or SuiteSparseQR_factorize (including the singleton-row R and the
361
// multifrontal R). A is m-by-n. B is n-by-nrhs, and X is m-by-nrhs.
362
363
template <typename Entry> static void rtsolve
364
(
365
// inputs
366
SuiteSparseQR_factorization <Entry> *QR,
367
int use_Q1fill,
368
369
Int nrhs, // number of columns of B
370
Int ldb, // leading dimension of B
371
Entry *B, // size n-by-nrhs with leading dimension ldb
372
373
// output, contents undefined on input
374
Entry *X, // size m-by-nrhs with leading dimension m
375
376
cholmod_common *cc
377
)
378
{
379
Entry xi ;
380
spqr_symbolic *QRsym ;
381
spqr_numeric <Entry> *QRnum ;
382
Int *R1p, *R1j, *Rmap, *Rp, *Rj, *Super, *HStair, *Hm, *Stair, *Q1fill,
383
*RmapInv ;
384
Entry *R1x, **Rblock, *R, *X1, *X2 ;
385
char *Rdead ;
386
Int i, j, k, m, n, p, kk, n1rows, n1cols, rank, nf, f, col1, col2, fp, pr,
387
fn, rm, row1, keepH, fm, h, t, live, jj ;
388
389
// -------------------------------------------------------------------------
390
// get the contents of the QR factorization
391
// -------------------------------------------------------------------------
392
393
QRsym = QR->QRsym ;
394
QRnum = QR->QRnum ;
395
n1rows = QR->n1rows ;
396
n1cols = QR->n1cols ;
397
n = QR->nacols ;
398
m = QR->narows ;
399
Q1fill = use_Q1fill ? QR->Q1fill : NULL ;
400
R1p = QR->R1p ;
401
R1j = QR->R1j ;
402
R1x = QR->R1x ;
403
Rmap = QR->Rmap ;
404
RmapInv = QR->RmapInv ;
405
rank = QR->rank ;
406
ASSERT ((Rmap == NULL) == (rank == n)) ;
407
ASSERT ((RmapInv == NULL) == (rank == n)) ;
408
ASSERT (rank <= MIN (m,n)) ;
409
410
keepH = QRnum->keepH ;
411
PR (("\nrtsolve keepH %ld n1rows %ld n1cols %ld\n", keepH, n1rows, n1cols));
412
413
// This code is only used when keepH is true. keepH is tested below,
414
// however, in case this changes in the future (the following assert would
415
// also need to be removed).
416
ASSERT (keepH) ;
417
418
nf = QRsym->nf ;
419
Rblock = QRnum->Rblock ;
420
Rp = QRsym->Rp ;
421
Rj = QRsym->Rj ;
422
Super = QRsym->Super ;
423
Rdead = QRnum->Rdead ;
424
425
HStair = QRnum->HStair ;
426
Hm = QRnum->Hm ;
427
428
ASSERT (ldb >= n) ;
429
430
// -------------------------------------------------------------------------
431
// X = E'*B or X = B
432
// -------------------------------------------------------------------------
433
434
X1 = X ;
435
if (rank == n)
436
{
437
// R is in upper triangular form
438
ASSERT (n <= m) ;
439
for (kk = 0 ; kk < nrhs ; kk++)
440
{
441
for (k = 0 ; k < n ; k++)
442
{
443
X1 [k] = B [Q1fill ? Q1fill [k] : k] ;
444
}
445
for ( ; k < m ; k++)
446
{
447
X1 [k] = 0 ;
448
}
449
X1 += m ;
450
B += ldb ;
451
}
452
}
453
else
454
{
455
// R is in squeezed form; use the mapping to trapezoidal
456
for (kk = 0 ; kk < nrhs ; kk++)
457
{
458
for (i = 0 ; i < rank ; i++)
459
{
460
k = RmapInv [i] ;
461
ASSERT (k >= 0 && k < n) ;
462
Int knew = Q1fill ? Q1fill [k] : k ;
463
ASSERT (knew >= 0 && knew < n) ;
464
X1 [i] = B [knew] ;
465
}
466
for ( ; i < m ; i++)
467
{
468
X1 [i] = 0 ;
469
}
470
X1 += m ;
471
B += ldb ;
472
}
473
}
474
475
// =========================================================================
476
// === solve with the singleton rows of R ==================================
477
// =========================================================================
478
479
X1 = X ;
480
for (kk = 0 ; kk < nrhs ; kk++)
481
{
482
for (i = 0 ; i < n1rows ; i++)
483
{
484
// divide by the diagonal (the singleton entry itself)
485
p = R1p [i] ;
486
xi = spqr_divide (X1 [i], spqr_conj (R1x [p]), cc) ;
487
X1 [i] = xi ;
488
489
// solve with the off-diagonal entries
490
for (p++ ; p < R1p [i+1] ; p++)
491
{
492
j = R1j [p] ;
493
ASSERT (j >= i && j < n) ;
494
k = Rmap ? Rmap [j] : j ;
495
ASSERT (k >= 0 && k < n) ;
496
if (k < rank)
497
{
498
// the jth column of R is live, and is the kth column of
499
// the squeezed R
500
X1 [k] -= spqr_conj (R1x [p]) * xi ;
501
}
502
}
503
}
504
X1 += m ;
505
}
506
507
// =========================================================================
508
// === solve with the multifrontal rows of R ===============================
509
// =========================================================================
510
511
Stair = NULL ;
512
fm = 0 ;
513
h = 0 ;
514
t = 0 ;
515
516
// start with row1 = n1rows, the first non-singleton row, which is the
517
// first row of the multifrontal R.
518
519
row1 = n1rows ;
520
for (f = 0 ; f < nf ; f++)
521
{
522
523
// ---------------------------------------------------------------------
524
// get the R block for front F
525
// ---------------------------------------------------------------------
526
527
R = Rblock [f] ;
528
col1 = Super [f] ; // first pivot column in front F
529
col2 = Super [f+1] ; // col2-1 is last pivot col
530
fp = col2 - col1 ; // number of pivots in front F
531
pr = Rp [f] ; // pointer to row indices for F
532
fn = Rp [f+1] - pr ; // # of columns in front F
533
534
if (keepH)
535
{
536
Stair = HStair + pr ; // staircase of front F
537
fm = Hm [f] ; // # of rows in front F
538
h = 0 ; // H vector starts in row h
539
}
540
541
// ---------------------------------------------------------------------
542
// solve with the squeezed upper triangular part of R
543
// ---------------------------------------------------------------------
544
545
rm = 0 ; // number of rows in R block
546
for (k = 0 ; k < fp ; k++)
547
{
548
j = col1 + k ;
549
ASSERT (Rj [pr + k] == j) ;
550
if (j+n1cols >= n) return ; // in case [A Binput] was factorized
551
if (keepH)
552
{
553
t = Stair [k] ; // length of R+H vector
554
ASSERT (t >= 0 && t <= fm) ;
555
if (t == 0)
556
{
557
live = FALSE ; // column k is dead
558
t = rm ; // dead col, R only, no H
559
h = rm ;
560
}
561
else
562
{
563
live = (rm < fm) ; // k is live, unless we hit the wall
564
h = rm + 1 ; // H vector starts in row h
565
}
566
ASSERT (t >= h) ;
567
}
568
else
569
{
570
// Ironically, this statement is itself dead code. In the
571
// current version, keepH is always TRUE whenever rtsolve is
572
// used. If H was not kept, this statement would determine
573
// whether the the pivot column is live or dead.
574
live = (!Rdead [j]) ; // DEAD
575
}
576
577
if (live)
578
{
579
// This is a live pivot column of the R block. rm has not been
580
// incremented yet, so at this point it is equal to the number
581
// of entries in the pivot column above the diagonal entry.
582
// The corresponding rows of the global R (and the solution X)
583
// are row1:row1+rm, and row1+rm is the "diagonal".
584
X1 = X + row1 ;
585
ASSERT (IMPLIES (Rmap != NULL, Rmap [j+n1cols] == row1 + rm)) ;
586
ASSERT (IMPLIES (Rmap == NULL, j+n1cols == row1 + rm)) ;
587
for (kk = 0 ; kk < nrhs ; kk++)
588
{
589
xi = X1 [rm] ; // xi = X (row1+rm,kk)
590
for (i = 0 ; i < rm ; i++)
591
{
592
xi -= spqr_conj (R [i]) * X1 [i] ;
593
}
594
X1 [rm] = spqr_divide (xi, spqr_conj (R [rm]), cc) ;
595
X1 += m ;
596
}
597
// After rm is incremented, it represents the number of rows
598
// so far in the R block
599
rm++ ;
600
}
601
602
// advance to the next column of R in the R block
603
R += rm + (keepH ? (t-h) : 0) ;
604
}
605
606
// ---------------------------------------------------------------------
607
// solve with the rectangular part of R
608
// ---------------------------------------------------------------------
609
610
for ( ; k < fn ; k++)
611
{
612
j = Rj [pr + k] ;
613
ASSERT (j >= col2 && j < QRsym->n) ;
614
jj = j + n1cols ;
615
if (jj >= n) break ; // in case [A Binput] was factorized
616
617
Int ii = Rmap ? Rmap [jj] : jj ;
618
ASSERT (ii >= n1rows && ii < n) ;
619
if (ii < rank)
620
{
621
// X (ii,:) -= R (row1:row+rm-1,:)' * X (row1:row1+rm-1,:)
622
X2 = X + row1 ; // pointer to X (row1,0)
623
X1 = X ;
624
for (kk = 0 ; kk < nrhs ; kk++)
625
{
626
xi = X1 [ii] ; // X (ii,kk)
627
for (i = 0 ; i < rm ; i++)
628
{
629
xi -= spqr_conj (R [i]) * X2 [i] ;
630
}
631
X1 [ii] = xi ;
632
X1 += m ;
633
X2 += m ;
634
}
635
}
636
637
// go to the next column of R
638
R += rm ;
639
if (keepH)
640
{
641
t = Stair [k] ; // length of R+H vector
642
ASSERT (t >= 0 && t <= fm) ;
643
h = MIN (h+1, fm) ; // H vector starts in row h
644
ASSERT (t >= h) ;
645
R += (t-h) ;
646
}
647
}
648
649
row1 += rm ; // advance past the rm rows of R in this front
650
}
651
ASSERT (row1 == rank) ;
652
}
653
654
655
// =============================================================================
656
// === SuiteSparseQR_solve (dense case) ========================================
657
// =============================================================================
658
659
// Use the R factor from a QR object returned by SuiteSparseQR_factorize to
660
// solve an upper or lower triangular system. SuiteSparseQR_factorize computes
661
// a QR object whose MATLAB equivalent is [Q,R,E]=spqr(A,opts) where opts.Q =
662
// 'Householder' (QR contains Q in Householder form, R, and the permutation
663
// vector E of 0:n-1). This function performs one of the following operations,
664
// where A is m-by-n:
665
//
666
// system=SPQR_RX_EQUALS_B (0): X = R\B B is m-by-k and X is n-by-k
667
// system=SPQR_RETX_EQUALS_B (1): X = E*(R\B) as above, E is a permutation
668
// system=SPQR_RTX_EQUALS_B (2): X = R'\B B is n-by-k and X is m-by-k
669
// system=SPQR_RTX_EQUALS_ETB (3): X = R'\(E'*B) as above, E is a permutation
670
//
671
// Both X and B are dense matrices. Only the first r rows and columns of R are
672
// used, where r is the rank estimate of A found by SuiteSparseQR_factorize.
673
674
template <typename Entry> cholmod_dense *SuiteSparseQR_solve // returns X
675
(
676
// inputs, not modified:
677
int system, // which system to solve (0,1,2,3)
678
SuiteSparseQR_factorization <Entry> *QR, // of an m-by-n sparse matrix A
679
cholmod_dense *B, // right-hand-side, m-by-nrhs or n-by-nrhs
680
// workspace and parameters
681
cholmod_common *cc
682
)
683
{
684
Entry *Bx ;
685
cholmod_dense *W, *X ;
686
Int m, n, nrhs, ldb, ok ;
687
688
// -------------------------------------------------------------------------
689
// get inputs
690
// -------------------------------------------------------------------------
691
692
RETURN_IF_NULL_COMMON (NULL) ;
693
RETURN_IF_NULL (B, NULL) ;
694
Int xtype = spqr_type <Entry> ( ) ;
695
RETURN_IF_XTYPE_INVALID (B, NULL) ;
696
RETURN_IF_NULL (QR, NULL) ;
697
RETURN_IF_NULL (QR->QRnum, NULL) ;
698
if (system < SPQR_RX_EQUALS_B || system > SPQR_RTX_EQUALS_ETB)
699
{
700
ERROR (CHOLMOD_INVALID, "Invalid system") ;
701
return (NULL) ;
702
}
703
m = QR->narows ;
704
n = QR->nacols ;
705
if ((Int) B->nrow != ((system <= SPQR_RETX_EQUALS_B) ? m : n))
706
{
707
ERROR (CHOLMOD_INVALID, "invalid dimensions") ;
708
return (NULL) ;
709
}
710
711
ASSERT (QR->bncols == 0) ;
712
713
cc->status = CHOLMOD_OK ;
714
715
nrhs = B->ncol ;
716
Bx = (Entry *) B->x ;
717
ldb = B->d ;
718
719
if (system == SPQR_RX_EQUALS_B || system == SPQR_RETX_EQUALS_B)
720
{
721
722
// ---------------------------------------------------------------------
723
// X = E*(R\B) or X=R\B
724
// ---------------------------------------------------------------------
725
726
Int *Rlive ;
727
Entry **Rcolp ;
728
X = cholmod_l_allocate_dense (n, nrhs, n, xtype, cc) ;
729
Int maxfrank = QR->QRnum->maxfrank ;
730
W = cholmod_l_allocate_dense (maxfrank, nrhs, maxfrank, xtype, cc) ;
731
Rlive = (Int *) cholmod_l_malloc (maxfrank, sizeof (Int), cc) ;
732
Rcolp = (Entry **) cholmod_l_malloc (maxfrank, sizeof (Entry *), cc) ;
733
ok = (X != NULL) && (W != NULL) && (cc->status == CHOLMOD_OK) ;
734
if (ok)
735
{
736
spqr_rsolve (QR, system == SPQR_RETX_EQUALS_B, nrhs, ldb, Bx,
737
(Entry *) X->x, Rcolp, Rlive, (Entry *) W->x, cc) ;
738
}
739
cholmod_l_free (maxfrank, sizeof (Int), Rlive, cc) ;
740
cholmod_l_free (maxfrank, sizeof (Entry *), Rcolp, cc) ;
741
cholmod_l_free_dense (&W, cc) ;
742
743
}
744
else
745
{
746
747
// ---------------------------------------------------------------------
748
// X = R'\(E'*B) or R'\B
749
// ---------------------------------------------------------------------
750
751
X = cholmod_l_allocate_dense (m, nrhs, m, xtype, cc) ;
752
ok = (X != NULL) ;
753
if (ok)
754
{
755
rtsolve (QR, system == SPQR_RTX_EQUALS_ETB, nrhs, ldb, Bx,
756
(Entry *) X->x, cc) ;
757
}
758
}
759
760
if (!ok)
761
{
762
// out of memory
763
ERROR (CHOLMOD_OUT_OF_MEMORY, "out of memory") ;
764
cholmod_l_free_dense (&X, cc) ;
765
return (NULL) ;
766
}
767
768
return (X) ;
769
}
770
771
772
template cholmod_dense *SuiteSparseQR_solve <Complex>
773
(
774
// inputs, not modified:
775
int system, // which system to solve
776
SuiteSparseQR_factorization <Complex> *QR, // of an m-by-n sparse matrix A
777
cholmod_dense *B, // right-hand-side, m-by-nrhs or n-by-nrhs
778
// workspace and parameters
779
cholmod_common *cc
780
) ;
781
782
template cholmod_dense *SuiteSparseQR_solve <double>
783
(
784
// inputs, not modified:
785
int system, // which system to solve
786
SuiteSparseQR_factorization <double> *QR, // of an m-by-n sparse matrix A
787
cholmod_dense *B, // right-hand-side, m-by-nrhs or n-by-nrhs
788
// workspace and parameters
789
cholmod_common *cc
790
) ;
791
792
793
// =============================================================================
794
// === SuiteSparseQR_solve (sparse case) =======================================
795
// =============================================================================
796
797
template <typename Entry> cholmod_sparse *SuiteSparseQR_solve // returns X
798
(
799
// inputs, not modified:
800
int system, // which system to solve (0,1,2,3)
801
SuiteSparseQR_factorization <Entry> *QR, // of an m-by-n sparse matrix A
802
cholmod_sparse *Bsparse, // right-hand-side, m-by-nrhs or n-by-nrhs
803
// workspace and parameters
804
cholmod_common *cc
805
)
806
{
807
cholmod_dense *Bdense, *Xdense ;
808
cholmod_sparse *Xsparse = NULL ;
809
RETURN_IF_NULL_COMMON (NULL) ;
810
RETURN_IF_NULL (QR, NULL) ;
811
RETURN_IF_NULL (Bsparse, NULL) ;
812
Int xtype = spqr_type <Entry> ( ) ;
813
RETURN_IF_XTYPE_INVALID (Bsparse, NULL) ;
814
cc->status = CHOLMOD_OK ;
815
816
Bdense = cholmod_l_sparse_to_dense (Bsparse, cc) ;
817
Xdense = SuiteSparseQR_solve <Entry> (system, QR, Bdense, cc) ;
818
cholmod_l_free_dense (&Bdense, cc) ;
819
Xsparse = cholmod_l_dense_to_sparse (Xdense, TRUE, cc) ;
820
cholmod_l_free_dense (&Xdense, cc) ;
821
822
if (Xsparse == NULL)
823
{
824
cc->status = CHOLMOD_OUT_OF_MEMORY ;
825
}
826
return (Xsparse) ;
827
}
828
829
template cholmod_sparse *SuiteSparseQR_solve <double>
830
(
831
// inputs, not modified:
832
int system, // which system to solve (0,1,2,3)
833
SuiteSparseQR_factorization <double> *QR, // of an m-by-n sparse matrix A
834
cholmod_sparse *Bsparse, // right-hand-side, m-by-nrhs or n-by-nrhs
835
// workspace and parameters
836
cholmod_common *cc
837
) ;
838
839
template cholmod_sparse *SuiteSparseQR_solve <Complex>
840
(
841
// inputs, not modified:
842
int system, // which system to solve (0,1,2,3)
843
SuiteSparseQR_factorization <Complex> *QR, // of an m-by-n sparse matrix A
844
cholmod_sparse *Bsparse, // right-hand-side, m-by-nrhs or n-by-nrhs
845
// workspace and parameters
846
cholmod_common *cc
847
) ;
848
849
850
// =============================================================================
851
// === get_H_vectors ===========================================================
852
// =============================================================================
853
854
// Get pointers to the Householder vectors in a single front.
855
// Returns # Householder vectors in F.
856
857
template <typename Entry> static Int get_H_vectors
858
(
859
// inputs
860
Int f, // front to operate on
861
SuiteSparseQR_factorization <Entry> *QR,
862
863
// outputs
864
Entry *H_Tau, // size QRsym->maxfn
865
Int *H_start, // size QRsym->maxfn
866
Int *H_end, // size QRsym->maxfn
867
cholmod_common *cc
868
)
869
{
870
spqr_symbolic *QRsym ;
871
spqr_numeric <Entry> *QRnum ;
872
Entry *Tau ;
873
Int *Rj, *Stair ;
874
Int col1, col2, fp, pr, fn, fm, h, nh, p, rm, k, j, t, n1cols, n ;
875
876
// -------------------------------------------------------------------------
877
// get the R block for front F
878
// -------------------------------------------------------------------------
879
880
QRsym = QR->QRsym ;
881
QRnum = QR->QRnum ;
882
n1cols = QR->n1cols ;
883
Rj = QRsym->Rj ;
884
n = QR->nacols ;
885
886
col1 = QRsym->Super [f] ; // first pivot column in front F
887
col2 = QRsym->Super [f+1] ; // col2-1 is last pivot col
888
fp = col2 - col1 ; // number of pivots in front F
889
pr = QRsym->Rp [f] ; // pointer to row indices for F
890
fn = QRsym->Rp [f+1] - pr ; // # of columns in front F
891
892
Stair = QRnum->HStair + pr ; // staircase of front F
893
Tau = QRnum->HTau + pr ; // Householder coeff. for front F
894
895
fm = QRnum->Hm [f] ; // # of rows in front F
896
h = 0 ; // H vector starts in row h
897
nh = 0 ; // number of Householder vectors
898
p = 0 ; // index into R+H block
899
900
// -------------------------------------------------------------------------
901
// find all the Householder vectors in the R+H block
902
// -------------------------------------------------------------------------
903
904
rm = 0 ; // number of rows in R block
905
for (k = 0 ; k < fn && nh < fm ; k++)
906
{
907
if (k < fp)
908
{
909
// a pivotal column of front F
910
j = col1 + k ;
911
ASSERT (Rj [pr + k] == j) ;
912
t = Stair [k] ; // length of R+H vector
913
ASSERT (t >= 0 && t <= fm) ;
914
if (t == 0)
915
{
916
p += rm ; // no H vector, skip over R
917
continue ;
918
}
919
else if (rm < fm)
920
{
921
rm++ ; // column k is not dead
922
}
923
h = rm ; // H vector starts in row h
924
ASSERT (t >= h) ;
925
}
926
else
927
{
928
// a non-pivotal column of front F
929
j = Rj [pr + k] ;
930
ASSERT (j >= col2 && j < QRsym->n) ;
931
t = Stair [k] ; // length of R+H vector
932
ASSERT (t >= rm && t <= fm) ;
933
h = MIN (h+1,fm) ; // one more row of C to skip
934
ASSERT (t >= h) ;
935
}
936
if (j+n1cols >= n) break ; // in case [A Binput] was factorized
937
p += rm ; // skip over R
938
H_Tau [nh] = Tau [k] ; // one more H vector found
939
H_start [nh] = p ; // H vector starts here
940
p += (t-h) ; // skip over the H vector
941
H_end [nh] = p ; // H vector ends here
942
nh++ ;
943
if (h == fm) break ; // that was the last H vector
944
}
945
946
return (nh) ;
947
}
948
949
950
// =============================================================================
951
// === load_H_vectors ==========================================================
952
// =============================================================================
953
954
// Load Householder vectors h1:h2-1 into the panel V. Return # of rows in V.
955
956
template <typename Entry> static Int load_H_vectors
957
(
958
// input
959
Int h1, // load vectors h1 to h2-1
960
Int h2,
961
Int *H_start, // vector h starts at R [H_start [h]]
962
Int *H_end, // vector h ends at R [H_end [h]-1]
963
Entry *R, // Rblock [f]
964
// output
965
Entry *V, // V is v-by-(h2-h1) and lower triangular
966
cholmod_common *cc
967
)
968
{
969
// v = length of last H vector
970
Int v = H_end [h2-1] - H_start [h2-1] + (h2-h1) ;
971
Entry *V1 = V ;
972
for (Int h = h1 ; h < h2 ; h++)
973
{
974
Int i ;
975
// This part of V is not accessed, for testing only:
976
// for (i = 0 ; i < h-h1 ; i++) V1 [i] = 0 ;
977
i = h-h1 ;
978
V1 [i++] = 1 ;
979
for (Int p = H_start [h] ; p < H_end [h] ; p++)
980
{
981
V1 [i++] = R [p] ;
982
}
983
for ( ; i < v ; i++)
984
{
985
V1 [i] = 0 ;
986
}
987
V1 += v ;
988
}
989
return (v) ;
990
}
991
992
993
// =============================================================================
994
// === Happly ==================================================================
995
// =============================================================================
996
997
// Given a QR factorization from spqr_1factor, apply the Householder vectors
998
// to a dense matrix X.
999
1000
template <typename Entry> static void Happly
1001
(
1002
// inputs
1003
int method, // 0,1,2,3
1004
SuiteSparseQR_factorization <Entry> *QR,
1005
Int hchunk, // apply hchunk Householder vectors at a time
1006
1007
// input/output
1008
Int m,
1009
Int n,
1010
Entry *X, // size m-by-n with leading dimension m; only
1011
// X (n1rows:m-1,:) or X (:,n1rows:n-1) is modified
1012
1013
// workspace, not defined on input or output
1014
Entry *H_Tau, // size QRsym->maxfn
1015
Int *H_start, // size QRsym->maxfn
1016
Int *H_end, // size QRsym->maxfn
1017
Entry *V, // size v-by-hchunk, where v = QRnum->maxfm
1018
Entry *C, // size: method 0,1: v*n, method 2,3: m*v
1019
Entry *W, // size: method 0,1: h*h+n*h, method 2,3: h*h+m*h
1020
// where h = hchunk
1021
cholmod_common *cc
1022
)
1023
{
1024
1025
spqr_symbolic *QRsym ;
1026
spqr_numeric <Entry> *QRnum ;
1027
Entry **Rblock, *R, *X2 ;
1028
Int *Hii, *Hip, *Hi ;
1029
Int nf, f, nh, h1, h2, v, n1rows, m2, n2 ;
1030
1031
// -------------------------------------------------------------------------
1032
// get the contents of the QR factorization
1033
// -------------------------------------------------------------------------
1034
1035
QRsym = QR->QRsym ;
1036
QRnum = QR->QRnum ;
1037
ASSERT (QRnum->keepH) ;
1038
nf = QRsym->nf ;
1039
Rblock = QRnum->Rblock ;
1040
Hii = QRnum->Hii ;
1041
Hip = QRsym->Hip ;
1042
ASSERT (QR->narows == ((method <= SPQR_QX) ? m : n)) ;
1043
n1rows = QR->n1rows ;
1044
1045
// -------------------------------------------------------------------------
1046
// operate on X (n1rows:m-1,:) or X (:,n1rows:n-1)
1047
// -------------------------------------------------------------------------
1048
1049
// The Householder vectors must be shifted downwards by n1rows, past the
1050
// singleton rows which are not part of the multifrontal part of the
1051
// QR factorization.
1052
1053
if (method == SPQR_QTX || method == SPQR_QX)
1054
{
1055
// Q*X or Q'*X
1056
X2 = X + n1rows ;
1057
n2 = n ;
1058
m2 = m - n1rows ;
1059
}
1060
else
1061
{
1062
X2 = X + n1rows * m ;
1063
n2 = n - n1rows ;
1064
m2 = m ;
1065
}
1066
1067
// -------------------------------------------------------------------------
1068
// apply the Householder vectors
1069
// -------------------------------------------------------------------------
1070
1071
if (method == SPQR_QTX || method == SPQR_XQ)
1072
{
1073
1074
// ---------------------------------------------------------------------
1075
// apply in forward direction
1076
// ---------------------------------------------------------------------
1077
1078
for (f = 0 ; f < nf ; f++)
1079
{
1080
// get the Householder vectors for front F
1081
nh = get_H_vectors (f, QR, H_Tau, H_start, H_end, cc) ;
1082
R = Rblock [f] ;
1083
Hi = &Hii [Hip [f]] ; // list of row indices of H
1084
1085
// apply the Householder vectors, one panel at a time
1086
for (h1 = 0 ; h1 < nh ; h1 = h2)
1087
{
1088
// load vectors h1:h2-1 from R into the panel V and apply them
1089
h2 = MIN (h1 + hchunk, nh) ;
1090
v = load_H_vectors (h1, h2, H_start, H_end, R, V, cc) ;
1091
ASSERT (v+h1 <= QR->QRnum->Hm [f]) ;
1092
spqr_panel (method, m2, n2, v, h2-h1, Hi+h1, V, H_Tau+h1,
1093
m, X2, C, W, cc) ;
1094
}
1095
}
1096
}
1097
else
1098
{
1099
1100
// ---------------------------------------------------------------------
1101
// apply in backward direction
1102
// ---------------------------------------------------------------------
1103
1104
for (f = nf-1 ; f >= 0 ; f--)
1105
{
1106
1107
// get the Householder vectors for front F
1108
nh = get_H_vectors (f, QR, H_Tau, H_start, H_end, cc) ;
1109
R = Rblock [f] ;
1110
Hi = &Hii [Hip [f]] ; // list of row indices of H
1111
1112
// apply the Householder vectors, one panel at a time
1113
for (h2 = nh ; h2 > 0 ; h2 = h1)
1114
{
1115
// load vectors h1:h2-1 from R into the panel V and apply them
1116
h1 = MAX (h2 - hchunk, 0) ;
1117
v = load_H_vectors (h1, h2, H_start, H_end, R, V, cc) ;
1118
ASSERT (v+h1 <= QR->QRnum->Hm [f]) ;
1119
spqr_panel (method, m2, n2, v, h2-h1, Hi+h1, V, H_Tau+h1,
1120
m, X2, C, W, cc) ;
1121
}
1122
}
1123
}
1124
}
1125
1126
1127
// =============================================================================
1128
// === SuiteSparseQR_qmult (dense case) ========================================
1129
// =============================================================================
1130
1131
// Applies Q in Householder form (as stored in the QR factorization object
1132
// returned by SuiteSparseQR_factorize) to a dense matrix X.
1133
//
1134
// method SPQR_QTX (0): Y = Q'*X
1135
// method SPQR_QX (1): Y = Q*X
1136
// method SPQR_XQT (2): Y = X*Q'
1137
// method SPQR_XQ (3): Y = X*Q
1138
1139
#define HCHUNK 32 // FUTURE: make this an input parameter
1140
1141
#define FREE_WORK \
1142
{ \
1143
cholmod_l_free_dense (&Zdense, cc) ; \
1144
cholmod_l_free_dense (&Vdense, cc) ; \
1145
cholmod_l_free_dense (&Wdense, cc) ; \
1146
cholmod_l_free_dense (&Cdense, cc) ; \
1147
cholmod_l_free (maxfn, sizeof (Entry), H_Tau, cc) ; \
1148
cholmod_l_free (maxfn, sizeof (Int), H_start, cc) ; \
1149
cholmod_l_free (maxfn, sizeof (Int), H_end, cc) ; \
1150
}
1151
1152
// returns Y of size m-by-n, or NULL on failure
1153
template <typename Entry> cholmod_dense *SuiteSparseQR_qmult
1154
(
1155
// inputs, not modified
1156
int method, // 0,1,2,3
1157
SuiteSparseQR_factorization <Entry> *QR,
1158
cholmod_dense *Xdense, // size m-by-n with leading dimension ldx
1159
1160
// workspace and parameters
1161
cholmod_common *cc
1162
)
1163
{
1164
cholmod_dense *Ydense, *Cdense, *Vdense, *Wdense, *Zdense ;
1165
Entry *X, *Y, *X1, *Y1, *Z1, *C, *V, *Z, *W, *H_Tau ;
1166
Int *HPinv, *H_start, *H_end ;
1167
Int i, k, mh, v, hchunk, ldx, m, n, maxfn, ok ;
1168
1169
// -------------------------------------------------------------------------
1170
// get inputs
1171
// -------------------------------------------------------------------------
1172
1173
RETURN_IF_NULL_COMMON (NULL) ;
1174
RETURN_IF_NULL (QR, NULL) ;
1175
RETURN_IF_NULL (QR->QRnum, NULL) ;
1176
RETURN_IF_NULL (QR->QRnum->Hm, NULL) ;
1177
RETURN_IF_NULL (Xdense, NULL) ;
1178
Int xtype = spqr_type <Entry> ( ) ;
1179
RETURN_IF_XTYPE_INVALID (Xdense, NULL) ;
1180
cc->status = CHOLMOD_OK ;
1181
1182
// get HPinv from QR->HP1inv if singletons exist, else QR->QRnum->HPinv
1183
HPinv = (QR->n1cols > 0) ? QR->HP1inv : QR->QRnum->HPinv ;
1184
1185
v = QR->QRnum->maxfm ;
1186
mh = QR->narows ;
1187
maxfn = QR->QRsym->maxfn ;
1188
1189
X = (Entry *) Xdense->x ;
1190
m = Xdense->nrow ;
1191
n = Xdense->ncol ;
1192
ldx = Xdense->d ;
1193
ASSERT (ldx >= m) ;
1194
1195
if (method == SPQR_QTX || method == SPQR_QX)
1196
{
1197
// rows of H and X must be the same
1198
if (mh != m)
1199
{
1200
ERROR (CHOLMOD_INVALID, "mismatched dimensions") ;
1201
return (NULL) ;
1202
}
1203
}
1204
else if (method == SPQR_XQT || method == SPQR_XQ)
1205
{
1206
// rows of H and columns of X must be the same
1207
if (mh != n)
1208
{
1209
ERROR (CHOLMOD_INVALID, "mismatched dimensions") ;
1210
return (NULL) ;
1211
}
1212
}
1213
else
1214
{
1215
ERROR (CHOLMOD_INVALID, "invalid method") ;
1216
return (NULL) ;
1217
}
1218
1219
// -------------------------------------------------------------------------
1220
// allocate result Y
1221
// -------------------------------------------------------------------------
1222
1223
Ydense = cholmod_l_allocate_dense (m, n, m, xtype, cc) ;
1224
if (cc->status < CHOLMOD_OK)
1225
{
1226
// out of memory
1227
return (NULL) ;
1228
}
1229
Y = (Entry *) Ydense->x ;
1230
1231
if (m == 0 || n == 0)
1232
{
1233
// nothing to do
1234
return (Ydense) ;
1235
}
1236
1237
// -------------------------------------------------------------------------
1238
// allocate workspace
1239
// -------------------------------------------------------------------------
1240
1241
Z = NULL ;
1242
Zdense = NULL ;
1243
ok = TRUE ;
1244
if (method == SPQR_QX || method == SPQR_XQT)
1245
{
1246
// Z of size m-by-n is needed only for Q*X and X*Q'
1247
Zdense = cholmod_l_allocate_dense (m, n, m, xtype, cc) ;
1248
ok = (Zdense != NULL) ;
1249
}
1250
1251
hchunk = HCHUNK ;
1252
ASSERT (v <= mh) ;
1253
1254
// C is workspace of size v-by-n or m-by-v
1255
Cdense = cholmod_l_allocate_dense (v, (method <= SPQR_QX) ? n : m,
1256
v, xtype, cc) ;
1257
Vdense = NULL ;
1258
Wdense = NULL ;
1259
1260
H_Tau = (Entry *) cholmod_l_malloc (maxfn, sizeof (Entry), cc) ;
1261
H_start = (Int *) cholmod_l_malloc (maxfn, sizeof (Int), cc) ;
1262
H_end = (Int *) cholmod_l_malloc (maxfn, sizeof (Int), cc) ;
1263
1264
if (!ok || Cdense == NULL || cc->status < CHOLMOD_OK)
1265
{
1266
// out of memory; free workspace and result Y
1267
ERROR (CHOLMOD_OUT_OF_MEMORY, "out of memory") ;
1268
cholmod_l_free_dense (&Ydense, cc) ;
1269
FREE_WORK ;
1270
return (NULL) ;
1271
}
1272
1273
// -------------------------------------------------------------------------
1274
// copy X into Z
1275
// -------------------------------------------------------------------------
1276
1277
if (method == SPQR_QX || method == SPQR_XQT)
1278
{
1279
// Z = X
1280
Z = (Entry *) Zdense->x ;
1281
Z1 = Z ;
1282
X1 = X ;
1283
for (k = 0 ; k < n ; k++)
1284
{
1285
for (i = 0 ; i < m ; i++)
1286
{
1287
Z1 [i] = X1 [i] ;
1288
}
1289
X1 += ldx ;
1290
Z1 += m ;
1291
}
1292
}
1293
1294
// -------------------------------------------------------------------------
1295
// allocate O(hchunk) workspace
1296
// -------------------------------------------------------------------------
1297
1298
// V is workspace of size v-by-hchunk
1299
Vdense = cholmod_l_allocate_dense (v, hchunk, v, xtype, cc) ;
1300
1301
// W is workspace of size h*h+n*h or h*h+m*h where h = hchunk
1302
Wdense = cholmod_l_allocate_dense (hchunk,
1303
hchunk + ((method <= SPQR_QX) ? n : m), hchunk, xtype, cc) ;
1304
1305
// -------------------------------------------------------------------------
1306
// punt if out of memory
1307
// -------------------------------------------------------------------------
1308
1309
if (Vdense == NULL || Wdense == NULL)
1310
{
1311
// PUNT: out of memory; try again with hchunk = 1
1312
cc->status = CHOLMOD_OK ;
1313
hchunk = 1 ;
1314
1315
cholmod_l_free_dense (&Vdense, cc) ;
1316
cholmod_l_free_dense (&Wdense, cc) ;
1317
1318
Vdense = cholmod_l_allocate_dense (v, hchunk, v, xtype, cc) ;
1319
Wdense = cholmod_l_allocate_dense (hchunk,
1320
hchunk + ((method <= SPQR_QX) ? n : m), hchunk, xtype, cc) ;
1321
1322
if (Vdense == NULL || Wdense == NULL)
1323
{
1324
// out of memory; free workspace and result Y
1325
ERROR (CHOLMOD_OUT_OF_MEMORY, "out of memory") ;
1326
cholmod_l_free_dense (&Ydense, cc) ;
1327
FREE_WORK ;
1328
return (NULL) ;
1329
}
1330
}
1331
1332
// the dimension (->nrow, ->ncol, and ->d) of this workspace is not used,
1333
// just the arrays themselves.
1334
V = (Entry *) Vdense->x ;
1335
C = (Entry *) Cdense->x ;
1336
W = (Entry *) Wdense->x ;
1337
1338
// -------------------------------------------------------------------------
1339
// Y = Q'*X, Q*X, X*Q, or X*Q'
1340
// -------------------------------------------------------------------------
1341
1342
PR (("Qfmult Method %ld m %ld n %ld X %p Y %p P %p\n", method, m, n, X, Y));
1343
1344
if (method == SPQR_QTX)
1345
{
1346
1347
// ---------------------------------------------------------------------
1348
// Y = Q'*X
1349
// ---------------------------------------------------------------------
1350
1351
// Y (P,:) = X, and change leading dimension from ldx to m
1352
X1 = X ;
1353
Y1 = Y ;
1354
for (k = 0 ; k < n ; k++)
1355
{
1356
for (i = 0 ; i < m ; i++)
1357
{
1358
Y1 [HPinv [i]] = X1 [i] ;
1359
}
1360
X1 += ldx ;
1361
Y1 += m ;
1362
}
1363
1364
// apply H to Y
1365
Happly (method, QR, hchunk, m, n, Y, H_Tau, H_start, H_end, V, C, W,cc);
1366
1367
}
1368
else if (method == SPQR_QX)
1369
{
1370
1371
// ---------------------------------------------------------------------
1372
// Y = Q*X
1373
// ---------------------------------------------------------------------
1374
1375
// apply H to Z
1376
Happly (method, QR, hchunk, m, n, Z, H_Tau, H_start, H_end, V, C, W,cc);
1377
1378
// Y = Z (P,:)
1379
Z1 = Z ;
1380
Y1 = Y ;
1381
for (k = 0 ; k < n ; k++)
1382
{
1383
for (i = 0 ; i < m ; i++)
1384
{
1385
Y1 [i] = Z1 [HPinv [i]] ;
1386
}
1387
Z1 += m ;
1388
Y1 += m ;
1389
}
1390
1391
}
1392
else if (method == SPQR_XQT)
1393
{
1394
1395
// ---------------------------------------------------------------------
1396
// Y = X*Q'
1397
// ---------------------------------------------------------------------
1398
1399
// apply H to Z
1400
Happly (method, QR, hchunk, m, n, Z, H_Tau, H_start, H_end, V, C, W,cc);
1401
1402
// Y = Z (:,P)
1403
Y1 = Y ;
1404
for (k = 0 ; k < n ; k++)
1405
{
1406
Z1 = Z + HPinv [k] * m ; // m = leading dimension of Z
1407
for (i = 0 ; i < m ; i++)
1408
{
1409
Y1 [i] = Z1 [i] ;
1410
}
1411
Y1 += m ;
1412
}
1413
1414
}
1415
else if (method == SPQR_XQ)
1416
{
1417
1418
// ---------------------------------------------------------------------
1419
// Y = X*Q
1420
// ---------------------------------------------------------------------
1421
1422
// Y (:,P) = X and change leading dimension from ldx to m
1423
X1 = X ;
1424
for (k = 0 ; k < n ; k++)
1425
{
1426
Y1 = Y + HPinv [k] * m ; // m = leading dimension of Y
1427
for (i = 0 ; i < m ; i++)
1428
{
1429
Y1 [i] = X1 [i] ;
1430
}
1431
X1 += ldx ;
1432
}
1433
1434
// apply H to Y
1435
Happly (method, QR, hchunk, m, n, Y, H_Tau, H_start, H_end, V, C, W,cc);
1436
1437
}
1438
1439
// -------------------------------------------------------------------------
1440
// free workspace and return Y
1441
// -------------------------------------------------------------------------
1442
1443
FREE_WORK ;
1444
1445
if (CHECK_BLAS_INT && !cc->blas_ok)
1446
{
1447
ERROR (CHOLMOD_INVALID, "problem too large for the BLAS") ;
1448
cholmod_l_free_dense (&Ydense, cc) ;
1449
return (NULL) ;
1450
}
1451
1452
return (Ydense) ;
1453
}
1454
1455
1456
template cholmod_dense *SuiteSparseQR_qmult <double>
1457
(
1458
// inputs, not modified
1459
int method, // 0,1,2,3
1460
SuiteSparseQR_factorization <double> *QR,
1461
cholmod_dense *Xdense, // size m-by-n with leading dimension ldx
1462
1463
// workspace and parameters
1464
cholmod_common *cc
1465
) ;
1466
1467
1468
template cholmod_dense *SuiteSparseQR_qmult <Complex>
1469
(
1470
// inputs, not modified
1471
int method, // 0,1,2,3
1472
SuiteSparseQR_factorization <Complex> *QR,
1473
cholmod_dense *Xdense, // size m-by-n with leading dimension ldx
1474
1475
// workspace and parameters
1476
cholmod_common *cc
1477
) ;
1478
1479
1480
// =============================================================================
1481
// === SuiteSparseQR_qmult (sparse case) =======================================
1482
// =============================================================================
1483
1484
// returns Y of size m-by-n, or NULL on failure
1485
template <typename Entry> cholmod_sparse *SuiteSparseQR_qmult
1486
(
1487
// inputs, not modified
1488
int method, // 0,1,2,3
1489
SuiteSparseQR_factorization <Entry> *QR,
1490
cholmod_sparse *Xsparse, // size m-by-n
1491
// workspace and parameters
1492
cholmod_common *cc
1493
)
1494
{
1495
cholmod_dense *Xdense, *Ydense ;
1496
cholmod_sparse *Ysparse = NULL ;
1497
RETURN_IF_NULL_COMMON (NULL) ;
1498
RETURN_IF_NULL (QR, NULL) ;
1499
RETURN_IF_NULL (Xsparse, NULL) ;
1500
Int xtype = spqr_type <Entry> ( ) ;
1501
RETURN_IF_XTYPE_INVALID (Xsparse, NULL) ;
1502
cc->status = CHOLMOD_OK ;
1503
1504
Xdense = cholmod_l_sparse_to_dense (Xsparse, cc) ;
1505
Ydense = SuiteSparseQR_qmult <Entry> (method, QR, Xdense, cc) ;
1506
cholmod_l_free_dense (&Xdense, cc) ;
1507
Ysparse = cholmod_l_dense_to_sparse (Ydense, TRUE, cc) ;
1508
cholmod_l_free_dense (&Ydense, cc) ;
1509
1510
if (Ysparse == NULL)
1511
{
1512
cc->status = CHOLMOD_OUT_OF_MEMORY ;
1513
}
1514
return (Ysparse) ;
1515
}
1516
1517
template cholmod_sparse *SuiteSparseQR_qmult <double>
1518
(
1519
// inputs, not modified
1520
int method, // 0,1,2,3
1521
SuiteSparseQR_factorization <double> *QR,
1522
cholmod_sparse *Xsparse, // size m-by-n
1523
// workspace and parameters
1524
cholmod_common *cc
1525
) ;
1526
1527
template cholmod_sparse *SuiteSparseQR_qmult <Complex>
1528
(
1529
// inputs, not modified
1530
int method, // 0,1,2,3
1531
SuiteSparseQR_factorization <Complex> *QR,
1532
cholmod_sparse *Xsparse, // size m-by-n
1533
// workspace and parameters
1534
cholmod_common *cc
1535
) ;
1536
1537
1538
// =============================================================================
1539
// === SuiteSparseQR_free ======================================================
1540
// =============================================================================
1541
1542
// Free the QR object; this is just a user-callable wrapper for spqr_freefac.
1543
1544
template <typename Entry> int SuiteSparseQR_free
1545
(
1546
SuiteSparseQR_factorization <Entry> **QR,
1547
cholmod_common *cc
1548
)
1549
{
1550
RETURN_IF_NULL_COMMON (FALSE) ;
1551
spqr_freefac <Entry> (QR, cc) ;
1552
return (TRUE) ;
1553
}
1554
1555
template int SuiteSparseQR_free <double>
1556
(
1557
SuiteSparseQR_factorization <double> **QR,
1558
cholmod_common *cc
1559
) ;
1560
1561
template int SuiteSparseQR_free <Complex>
1562
(
1563
SuiteSparseQR_factorization <Complex> **QR,
1564
cholmod_common *cc
1565
) ;
1566
1567
1568
// =============================================================================
1569
// === SuiteSparseQR_min2norm ==================================================
1570
// =============================================================================
1571
1572
// Find the min 2-norm solution for underdetermined systems (A is m-by-n with
1573
// m<n), or find a least-squares solution otherwise.
1574
1575
template <typename Entry> cholmod_dense *SuiteSparseQR_min2norm
1576
(
1577
int ordering, // all, except 3:given treated as 0:fixed
1578
double tol,
1579
cholmod_sparse *A,
1580
cholmod_dense *B,
1581
cholmod_common *cc
1582
)
1583
{
1584
cholmod_dense *X ;
1585
RETURN_IF_NULL_COMMON (NULL) ;
1586
RETURN_IF_NULL (A, NULL) ;
1587
RETURN_IF_NULL (B, NULL) ;
1588
Int xtype = spqr_type <Entry> ( ) ;
1589
RETURN_IF_XTYPE_INVALID (A, NULL) ;
1590
RETURN_IF_XTYPE_INVALID (B, NULL) ;
1591
cc->status = CHOLMOD_OK ;
1592
1593
if (A->nrow < A->ncol)
1594
{
1595
#ifdef TIMING
1596
double t0 = spqr_time ( ) ;
1597
#endif
1598
1599
// x=A\B, using a QR factorization of A'
1600
SuiteSparseQR_factorization <Entry> *QR ;
1601
cholmod_sparse *AT ;
1602
cholmod_dense *Y ;
1603
// [Q,R,E] = qr (A')
1604
AT = cholmod_l_transpose (A, 2, cc) ;
1605
QR = SuiteSparseQR_factorize <Entry> (ordering, tol, AT, cc);
1606
cholmod_l_free_sparse (&AT, cc) ;
1607
// solve Y = R'\(E'*B)
1608
Y = SuiteSparseQR_solve (SPQR_RTX_EQUALS_ETB, QR, B, cc) ;
1609
// X = Q*Y
1610
X = SuiteSparseQR_qmult (SPQR_QX, QR, Y, cc) ;
1611
cholmod_l_free_dense (&Y, cc) ;
1612
spqr_freefac (&QR, cc) ;
1613
1614
#ifdef TIMING
1615
double t3 = spqr_time ( ) ;
1616
double total_time = t3 - t0 ;
1617
cc->other1 [3] = total_time - cc->other1 [1] - cc->other1 [2] ;
1618
#endif
1619
1620
}
1621
else
1622
{
1623
// x=A\B, using a QR factorization of A
1624
SuiteSparseQR <Entry> (ordering, tol, 0, 2, A, NULL, B, NULL, &X,
1625
NULL, NULL, NULL, NULL, NULL, cc) ;
1626
}
1627
1628
if (X == NULL)
1629
{
1630
// After A and B are valid, the only way any of the above functions can
1631
// fail is by running out of memory. However, if (say) AT is NULL
1632
// because of insufficient memory, SuiteSparseQR_factorize will report
1633
// CHOLMOD_INVALID since its input is NULL. In that case, restore the
1634
// error back to CHOLMOD_OUT_OF_MEMORY.
1635
cc->status = CHOLMOD_OUT_OF_MEMORY ;
1636
}
1637
1638
return (X) ;
1639
}
1640
1641
template cholmod_dense *SuiteSparseQR_min2norm <double>
1642
(
1643
int ordering, // all, except 3:given treated as 0:fixed
1644
double tol,
1645
cholmod_sparse *A,
1646
cholmod_dense *B,
1647
cholmod_common *cc
1648
) ;
1649
1650
template cholmod_dense *SuiteSparseQR_min2norm <Complex>
1651
(
1652
int ordering, // all, except 3:given treated as 0:fixed
1653
double tol,
1654
cholmod_sparse *A,
1655
cholmod_dense *B,
1656
cholmod_common *cc
1657
) ;
1658
1659
1660
// =============================================================================
1661
// === SuiteSparseQR_min2norm (sparse case) ====================================
1662
// =============================================================================
1663
1664
template <typename Entry> cholmod_sparse *SuiteSparseQR_min2norm // returns X
1665
(
1666
int ordering, // all, except 3:given treated as 0:fixed
1667
double tol,
1668
cholmod_sparse *A,
1669
cholmod_sparse *Bsparse,
1670
cholmod_common *cc
1671
)
1672
{
1673
#ifdef TIMING
1674
double t0 = spqr_time ( ) ;
1675
#endif
1676
1677
cholmod_dense *Bdense, *Xdense ;
1678
cholmod_sparse *Xsparse = NULL ;
1679
RETURN_IF_NULL_COMMON (NULL) ;
1680
RETURN_IF_NULL (A, NULL) ;
1681
RETURN_IF_NULL (Bsparse, NULL) ;
1682
Int xtype = spqr_type <Entry> ( ) ;
1683
RETURN_IF_XTYPE_INVALID (A, NULL) ;
1684
RETURN_IF_XTYPE_INVALID (Bsparse, NULL) ;
1685
cc->status = CHOLMOD_OK ;
1686
1687
Bdense = cholmod_l_sparse_to_dense (Bsparse, cc) ;
1688
Xdense = SuiteSparseQR_min2norm <Entry> (ordering, tol, A, Bdense, cc) ;
1689
cholmod_l_free_dense (&Bdense, cc) ;
1690
Xsparse = cholmod_l_dense_to_sparse (Xdense, TRUE, cc) ;
1691
cholmod_l_free_dense (&Xdense, cc) ;
1692
1693
if (Xsparse == NULL)
1694
{
1695
cc->status = CHOLMOD_OUT_OF_MEMORY ;
1696
}
1697
1698
#ifdef TIMING
1699
double t3 = spqr_time ( ) ;
1700
double total_time = t3 - t0 ;
1701
cc->other1 [3] = total_time - cc->other1 [1] - cc->other1 [2] ;
1702
#endif
1703
1704
return (Xsparse) ;
1705
}
1706
1707
template cholmod_sparse *SuiteSparseQR_min2norm <double>
1708
(
1709
int ordering, // all, except 3:given treated as 0:fixed
1710
double tol,
1711
cholmod_sparse *A,
1712
cholmod_sparse *Bsparse,
1713
cholmod_common *cc
1714
) ;
1715
1716
template cholmod_sparse *SuiteSparseQR_min2norm <Complex>
1717
(
1718
int ordering, // all, except 3:given treated as 0:fixed
1719
double tol,
1720
cholmod_sparse *A,
1721
cholmod_sparse *Bsparse,
1722
cholmod_common *cc
1723
) ;
1724
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1