#
source:
trunk/projects/slim-curve/src/main/c/GCI_Lsqnonneg.c
@
7607

Revision 7607, 19.5 KB checked in by paulbarber, 9 years ago (diff) |
---|

Line | |
---|---|

1 | // Updated file from JG, 2.11.02 |

2 | |

3 | /* Non-neg least squares library function |

4 | |

5 | This code is based on code and descriptions in the book |

6 | Lawson, C.L. and Hanson, R.J., Solving Least Squares Problems, |

7 | Prentice-Hall, 1974 |

8 | |

9 | A brief description of the various algorithms used is included |

10 | within this file, but for a full explanation, please refer to the |

11 | above-mentioned book or a similar reference. |

12 | |

13 | Julian Gilbey, Gray Cancer Institute, September 2002 |

14 | */ |

15 | |

16 | /* Include files we need */ |

17 | #include <stdio.h> /* for NULL */ |

18 | #include <stdlib.h> /* for malloc/free */ |

19 | #include <string.h> /* for memcpy */ |

20 | #include <math.h> /* for fabs and sqrt */ |

21 | #ifndef max |

22 | #define max(a,b) ((a)<(b) ? (b) : (a)) |

23 | #endif |

24 | |

25 | /* The greatest number of variables to be determined. This is to |

26 | allocate an array in advance to save having to use malloc/free */ |

27 | #define MAX_VARS 50 |

28 | |

29 | /* The greatest number of equations to be used. Again, this is |

30 | to allocate an array in advance */ |

31 | #define MAX_EQNS 10000 |

32 | |

33 | /* Some tolerence calculations */ |

34 | /* Approximate machine epsilon */ |

35 | #define EPS 1e-9 |

36 | #define TOL (100*EPS) |

37 | |

38 | |

39 | /********************************************************************** |

40 | Function Householder |

41 | |

42 | This function takes as input a vector, a pivot element and a zero |

43 | point, and applies a Householder transformation to the vector to |

44 | make all of its entries zero from the point specified, pivoting |

45 | around the specified pivot entry, modifying the vector in the |

46 | process, as described below. This function can also take a |

47 | collection of other vectors and apply this Householder |

48 | transformation to these vectors. |

49 | |

50 | Here is a brief description of Householder transformations: we are |

51 | looking for an orthogonal transformation Q with |

52 | |

53 | Qv = ( v_0 v_1 ... v_{p-1} -sigma.sqrt(v_p^2+sum_{i=l}^{m-1} v_i^2) |

54 | v_{p+1} ... v_{l-1} 0 ... 0)^T =: y |

55 | |

56 | where sigma = +1 if v_p >= 0 and -1 if v_p < 0. |

57 | |

58 | Here p is the pivot index and l is the zeroing index, where we |

59 | assume 0 <= p <= l < m (where v is an m-vector, zero-based). If |

60 | these assumptions are not true, we take Q=I. |

61 | |

62 | Q is given by Q = I - 2uu^T/|u|^2 for some vector u, but we do |

63 | not actually need to calculate it. We do the following: |

64 | |

65 | Compute: s := -sigma.sqrt(v_p^2+sum_{i=l}^m v_i^2) |

66 | u_i := 0, i=0, ..., p-1 |

67 | u_p := v_p - s |

68 | u_i := 0, i=p+1, ..., l-1 |

69 | u_i := v_i, i=l, ..., m-1 |

70 | b := su_p (so b = -|u|^2/2) |

71 | Q := ( I_m + b^{-1}uu^T if b != 0 |

72 | ( I_m if b = 0 |

73 | |

74 | After this, we place the non-zero components of u and y into the |

75 | storage previously occupied by v; we let v_p := y_p and store u_p |

76 | separately. |

77 | |

78 | To apply this transformation to the set of m-vectors c_j, getting |

79 | d_j = Qc_j, we do: |

80 | |

81 | t_j := b^{-1}(u^T c_j) |

82 | d_j := c_j + t_j u |

83 | |

84 | if b != 0, and d_j := c_j if b = 0. |

85 | |

86 | This function computes u, b, y:=Qv and d_j := Qc_j for all c_j we |

87 | are given. If u and y have already been computed, we can apply |

88 | this Q to new c_j vectors. |

89 | |

90 | Arguments: mode: 1 means calculate u, b, y and then Qc_j |

91 | 2 means u and y are precalculated, just calculate Qc_j |

92 | v, the vector to pivot |

93 | NB this vector is modified during execution of this |

94 | function, as described above |

95 | p, the pivot |

96 | l, the zeroing point |

97 | m, the length of the vector v |

98 | u_p (double *), the calculated value of u_p |

99 | C, double **C, array of pointers to the c_j vectors |

100 | nC, number of c_j vectors |

101 | If C=NULL or nC=0, the code will not attempt to find any |

102 | Qc_j. |

103 | NB the calculated d_j := Qc_j vectors will overwrite the |

104 | C array |

105 | |

106 | Here is the collected algorithm, as outlined above: |

107 | |

108 | 0 // Start here if mode = 1 |

109 | 1 s := sqrt(v_p^2 + sum_{i=l}^{m-1} v_i^2) |

110 | (Use the calculation trick described below under Givens |

111 | transformations to avoid calculation errors) |

112 | 2 if v_p > 0, s := -s |

113 | 3 u_p := v_p - s, v_p := s |

114 | 4 // Start here if mode = 2 |

115 | 5 b := v_p.u_p |

116 | 6 if b=0 or no C, stop |

117 | 7 for each c_j do steps 8-10 |

118 | 8 t_j := (c_j[p].u_p + sum_{i=l}^{m-1} c_j[i].v_i) / b |

119 | 9 c_j[p] := c_j[p] + t_j.u_p |

120 | 10 for i=l,...,m-1, c_j[i] := c_j[i] + t_j.v_i |

121 | |

122 | */ |

123 | |

124 | void Householder(int mode, double *v, int p, int l, int m, double *u_p, |

125 | double *C[], int nC) |

126 | { |

127 | double vmax, vmax_inv, vtmp, s, b, tj; |

128 | double *cj; |

129 | int i, j; |

130 | |

131 | if (0>p || p>=l || l>m) return; |

132 | |

133 | /* find the maximum v_i */ |

134 | vmax = fabs(v[p]); |

135 | if (mode == 1) { |

136 | for (i=l; i<m; i++) |

137 | vmax = max(vmax, fabs(v[i])); |

138 | if (vmax == 0) return; /* oops; not much to do here, really */ |

139 | vmax_inv = 1/vmax; /* only calculate this once; does it matter? */ |

140 | |

141 | /* we now calculate s, which is step 1 in the algorithm */ |

142 | vtmp = v[p]*vmax_inv; |

143 | s = vtmp*vtmp; |

144 | for (i=l; i<m; i++) { |

145 | vtmp = v[i]*vmax_inv; |

146 | s += vtmp*vtmp; |

147 | } |

148 | s = vmax*sqrt(s); |

149 | /* this is step 1 done */ |

150 | if (v[p]>0) s = -s; |

151 | *u_p = v[p]-s; |

152 | v[p] = s; |

153 | } else { |

154 | if (vmax == 0) return; /* same test as before, only now we are |

155 | testing only the v_p element, which is |

156 | actually y_p */ |

157 | } |

158 | |

159 | /* Now calculate the Qc_j vectors */ |

160 | if (nC <= 0 || C==NULL) return; |

161 | b = v[p] * (*u_p); |

162 | if (b >= 0) return; |

163 | |

164 | for (j=0; j<nC; j++) { |

165 | cj = C[j]; |

166 | tj = cj[p] * (*u_p); |

167 | for (i=l; i<m; i++) |

168 | tj += cj[i]*v[i]; |

169 | |

170 | if (tj != 0) { |

171 | tj /= b; |

172 | cj[p] += tj*(*u_p); |

173 | for (i=l; i<m; i++) |

174 | cj[i] += tj*v[i]; |

175 | } |

176 | } |

177 | } |

178 | |

179 | |

180 | /********************************************************************** |

181 | Functions GivensCalc, GivensApply |

182 | |

183 | The GivensCalc function takes as input a 2-vector, and computes |

184 | the Givens rotation for this vector. |

185 | |

186 | If the vector is v = (v_1 v_2)^T != 0, then the Givens rotation is |

187 | the 2x2 (special) orthogonal matrix |

188 | |

189 | G = ( c s ) |

190 | ( -s c ) |

191 | |

192 | (with c^2 + s^2 = 1) such that Gv = (r 0)^T, where |

193 | r=|v|=sqrt(v_1^2 + v_2^2). |

194 | |

195 | This function calculates c, s and r, avoiding overflow and |

196 | underflow issues, and returns these values. (We can allow r to |

197 | overwrite v_1 or v_2 if we wish.) The calculation is as follows: |

198 | |

199 | To compute r=sqrt(x^2 + y^2): |

200 | t := max {|x|, |y|} |

201 | u := min {|x|, |y|} |

202 | r := ( t.sqrt(1+(u/t)^2) if t != 0 |

203 | ( 0 if t = 0 |

204 | |

205 | The function GivensApply applies a Givens rotation to a |

206 | 2-vector. No rocket science there. |

207 | |

208 | Here, then, are the algorithms we use: |

209 | |

210 | Calculate a Givens rotation for the vector (v_1 v_2)^T |

211 | (GivensCalc): |

212 | |

213 | Arguments: v_1, v_2, the input vector |

214 | c, s, r, pointers to the answers. r may point to v_1 or v_2 |

215 | |

216 | 1 if |v_1| <= |v_2| go to step 8 |

217 | 2 w := v_2/v_1 |

218 | 3 q := sqrt(1+w^2) |

219 | 4 c := 1/q |

220 | 5 if v_1 < 0, c := -c |

221 | 6 s := wc |

222 | 7 r := |v_1|.q and stop |

223 | 8 if v_2 != 0, go to step 10 |

224 | 9 [v_1=v_2=0] c := 1, s:= 0, r := 0, stop |

225 | 10 w := v_1/v_2 |

226 | 11 q := sqrt(1+w^2) |

227 | 12 s := 1/q |

228 | 13 if v_1 < 0, s := -s |

229 | 14 c := ws |

230 | 15 r := |v_2|.q and stop |

231 | |

232 | Apply the Givens rotation G to the vector (z_1 z_2)^T |

233 | (GivensApply): |

234 | |

235 | Arguments: c, s, components of the Givens rotation matrix |

236 | z_1, z_2, (double *) the vector to apply it to; these |

237 | will be overwritten! |

238 | |

239 | 1 w := z_1.c + z_2.s |

240 | 2 z_2 := -z_1.s + z_2.c |

241 | 3 z_1 := w |

242 | |

243 | */ |

244 | |

245 | void GivensCalc(double v_1, double v_2, double *c, double *s, double *r) |

246 | { |

247 | double w, q; |

248 | |

249 | if (fabs(v_1) > fabs(v_2)) { |

250 | w = v_2/v_1; |

251 | q = sqrt(1+w*w); |

252 | *c = 1/q; |

253 | if (v_1 < 0) *c = - (*c); |

254 | *s = w*(*c); |

255 | *r = fabs(v_1)*q; |

256 | } else { |

257 | if (v_2 == 0) { |

258 | *c = 1; *s = 0; *r = 0; |

259 | } else { |

260 | w = v_1/v_2; |

261 | q = sqrt(1+w*w); |

262 | *s = 1/q; |

263 | if (v_2 < 0) *s = - (*s); |

264 | *c = w*(*s); |

265 | *r = fabs(v_2)*q; |

266 | } |

267 | } |

268 | } |

269 | |

270 | |

271 | void GivensApply(double c, double s, double *z_1, double *z_2) |

272 | { |

273 | double w; |

274 | |

275 | w = (*z_1)*c + (*z_2)*s; |

276 | *z_2 = -(*z_1)*s + (*z_2)*c; |

277 | *z_1 = w; |

278 | } |

279 | |

280 | |

281 | /********************************************************************** |

282 | Function GCI_lsqnonneg |

283 | |

284 | This function solves the non-negative least squares problem: |

285 | minimise |Ax-b| subject to x >= 0 (where |v| is the 2-norm of the |

286 | vector v). |

287 | |

288 | Arguments: A, an m x n matrix, in the form double A[n][m], so |

289 | the columns of A are A[0], A[1], ..., A[n-1] |

290 | b, the m-vector |

291 | |

292 | !!! NB: A and B will both be modified unless preserve is |

293 | !!! non-zero (see below). |

294 | |

295 | x, the solution will be placed here |

296 | m ) the dimensions of A |

297 | n ) |

298 | preserve, copy A and b before solving the problem |

299 | rnorm, (double *) the value of |Ax-b| with the |

300 | determined x if the function was successful or if the |

301 | iteration count was exceeded. This can be NULL. |

302 | lambda, an n-vector which will contain the dual vector |

303 | on completion (that is, the Lagrange multipliers). |

304 | This can be NULL. |

305 | |

306 | On exit: |

307 | The return value will be 0 on success, and negative if |

308 | a problem occurred: |

309 | -1: m > MAX_EQNS or m <= 0 |

310 | -2: n > MAX_VARS or n <= 0 |

311 | -3: iteration count exceeded: more than 3*n iterations |

312 | performed |

313 | -4: memory allocation problems |

314 | |

315 | The algorithm is similar to the simplex algorithm, and uses |

316 | Lagrange multipliers. |

317 | |

318 | 1 Set P:={}, Z:={1,2,...,n}, x:=0 |

319 | 2 Compute w := A^T (b-Ax) |

320 | 3 If Z={} or if w_j <= 0 for all j in Z, stop |

321 | 4 Find t in Z such that w_t = max { w_j : j in Z } |

322 | 5 Move index t from Z to P |

323 | 6 Let A_P denote the m x n matrix defined by |

324 | |

325 | ( column j of A if j in P |

326 | column j of A_P := ( |

327 | ( 0 if j in Z |

328 | |

329 | Compute the n-vector z as a solution of the regular least squares |

330 | problem minimise | A_P z - b |. Note that only the components |

331 | z_j, j in P, are determined by this problem. Define z_j:=0 for j in Z. |

332 | 7 If z_j>0 for all j in P, then set x:=z and goto step 2 |

333 | 8 Find an index q in P such that |

334 | x_q/(x_q-z_q) = min { x_j/(x_j-z_j) : z_j <= 0, j in P } |

335 | 9 Set alpha := x_q/(x_q-z_q) |

336 | 10 Set x := x + alpha(z-x) |

337 | 11 Move from set P to set Z all j in P for which x_j=0. Go to step 6. |

338 | |

339 | On termination, x satisfies x_j>0 for j in P and x_j=0 for j in Z, and |

340 | x is a solution to the least squares problem minimise |A_P z - b|. |

341 | |

342 | Step 6 (finding the solution to the least squares problem |

343 | |A_P z - b|) is performed using a QR factorisation of A_P, which in |

344 | turn is calculated using Householder matrices and Givens rotations, |

345 | and clever tricks for keeping track of the QR factorisation of A_P |

346 | as P changes, as we describe below. |

347 | */ |

348 | |

349 | #define Amatrix(row,col) A[col][row] |

350 | |

351 | int GCI_lsqnonneg(double **A_orig, double *b_orig, double *x, int m, int n, |

352 | int preserve, double *rnorm_orig, double *lambda) |

353 | { |

354 | double *AA[MAX_VARS], bb[MAX_EQNS], **A, *b; /* to preserve originals */ |

355 | double rnorm2, *rnorm; /* in case rnorm_orig is NULL */ |

356 | double ww[MAX_VARS], z[MAX_EQNS], *w; |

357 | double *C[MAX_VARS]; /* for passing to Householder */ |

358 | double wmax, asave, u_p, unorm, sum; /* miscellany */ |

359 | int iter, itmax; |

360 | int index[MAX_VARS], p; |

361 | int i, ii, j, jj, k, l, q, zmax, index_zmax; |

362 | |

363 | /* Meanings of variables: |

364 | index[] is an array containing information on the sets Z and P. |

365 | index[0..p-1] is the set P |

366 | index[p..n-1] is the set Z |

367 | |

368 | z[] is working space |

369 | */ |

370 | |

371 | /* Check variables */ |

372 | if (m > MAX_EQNS || m <= 0) return -1; |

373 | if (n > MAX_VARS || n <= 0) return -2; |

374 | w = (lambda == NULL) ? ww : lambda; |

375 | rnorm = (rnorm_orig == NULL) ? &rnorm2 : rnorm_orig; |

376 | if (preserve) { |

377 | /* We allocate one long array and split it into pieces */ |

378 | if ((AA[0] = malloc(m*n*sizeof(double))) == NULL) |

379 | return -4; |

380 | for (i=0; i<n; i++) { |

381 | AA[i] = AA[0]+m*i; |

382 | for (j=0; j<m; j++) |

383 | AA[i][j] = A_orig[i][j]; |

384 | } |

385 | for (j=0; j<m; j++) |

386 | bb[j] = b_orig[j]; |

387 | A = AA; |

388 | b = bb; |

389 | } else { |

390 | A = A_orig; |

391 | b = b_orig; |

392 | } |

393 | |

394 | iter = 0; |

395 | itmax = 3*n; |

396 | |

397 | /* Initialise the index and x arrays */ |

398 | p = 0; |

399 | for (i=0; i<n; i++) { |

400 | x[i] = 0; |

401 | index[i] = i; |

402 | } |

403 | |

404 | /* Main loop */ |

405 | while (p<n && p<m) { |

406 | /* Compute the dual (negative gradient) vector w=A^T(b-Ax) |

407 | We are only interested in the components of this vector in Z, |

408 | so we only calculate these. It turns out that for these |

409 | components, we have w_j=(A^T.b)_j, and that (A^T.Ax)_j=0, |

410 | remembering that A is continually modified. I really don't |

411 | yet know why this is true. */ |

412 | |

413 | for (i=p; i<n; i++) { |

414 | j = index[i]; |

415 | sum = 0; |

416 | for (k=p; k<m; k++) |

417 | sum += Amatrix(k,j) * b[k]; |

418 | w[j] = sum; |

419 | } |

420 | |

421 | wmax = 0; |

422 | while (wmax == 0) { |

423 | /* We now find the maximum positive w_j. */ |

424 | zmax = -1; |

425 | for (i=p; i<n; i++) { |

426 | j = index[i]; |

427 | if (w[j]>wmax) { |

428 | wmax = w[j]; |

429 | zmax = i; |

430 | } |

431 | } |

432 | |

433 | if (wmax <= 0) /* all w_j non-positive, so done */ |

434 | goto terminate; |

435 | |

436 | /* We move index[zmax] from Z to P |

437 | We begin the transformation and check the new diagonal element |

438 | to avoid near linear dependence */ |

439 | index_zmax = index[zmax]; |

440 | asave = Amatrix(p,index_zmax); |

441 | |

442 | Householder(1, A[index_zmax], p, p+1, m, &u_p, NULL, 0); |

443 | /* u has ended up in the index_zmax-th column of A */ |

444 | unorm = 0; |

445 | for (k=0; k<p; k++) |

446 | unorm += Amatrix(k,index_zmax)*Amatrix(k,index_zmax); |

447 | unorm = sqrt(unorm); |

448 | |

449 | /* Is the index_zmax-th column sufficiently independent |

450 | (whatever that means)? */ |

451 | if ((unorm != 0 && fabs(Amatrix(p,index_zmax)/unorm) < TOL) || |

452 | (unorm == 0 && Amatrix(p,index_zmax) == 0)) { |

453 | /* no, it isn't */ |

454 | Amatrix(p,index_zmax) = asave; |

455 | w[index_zmax] = 0; |

456 | wmax = 0; |

457 | } else { |

458 | /* it is: copy b[] into z[], update z[] and solve for ztest, |

459 | the proposed new value for x[index_zmax] */ |

460 | double ztest; |

461 | for (i=0; i<m; i++) |

462 | z[i] = b[i]; |

463 | C[0] = z; /* Need to pass Householder an array of pointers */ |

464 | Householder(2, A[index_zmax], p, p+1, m, &u_p, C, 1); |

465 | ztest = z[p] / Amatrix(p,index_zmax); |

466 | if (ztest <= 0) { |

467 | /* nah, this won't do, try a different zmax */ |

468 | Amatrix(p,index_zmax) = asave; |

469 | w[index_zmax] = 0; |

470 | wmax = 0; |

471 | } |

472 | /* We've now got an acceptable index to move from Z to P! |

473 | We'll leave this while (wmax==0) loop now, then. */ |

474 | } |

475 | } |

476 | |

477 | /* And take heed of the fact that we've got our index */ |

478 | for (i=0; i<m; i++) |

479 | b[i] = z[i]; |

480 | index[zmax] = index[p]; |

481 | index[p] = index_zmax; |

482 | p++; |

483 | |

484 | /* We now apply our Householder transformation to all of the |

485 | columns of A which correspond to index[p], ..., index[n-1]. We |

486 | could call the Householder function for each column separately, |

487 | but it's probably more efficient (fewer function calls) to |

488 | bundle the columns up into one new temporary array. Note that |

489 | we've incremented p! */ |

490 | if (p<n) { |

491 | for (i=p; i<n; i++) { |

492 | C[i-p] = A[index[i]]; /* the index[i]-th column of A */ |

493 | } |

494 | Householder(2, A[index_zmax], p-1, p, m, &u_p, C, n-p); |

495 | } |

496 | |

497 | /* ... zero the index_zmax-th column of A below the |

498 | diagonal, and set w[index_zmax]=0 ... */ |

499 | for (i=p; i<m; i++) |

500 | Amatrix(i,index_zmax) = 0; |

501 | w[index_zmax] = 0; |

502 | |

503 | /* ... and solve the triangular system to solve the least squares |

504 | problem |A_P.z-b| (step 6 of algorithm), putting the solution |

505 | into z. Note also that a copy of b is still stored in z when |

506 | we begin. |

507 | |

508 | We have the matrix A_P which has j-th column being index[j], |

509 | j=0, ..., p-1. The j-th column is zero below the diagonal. So |

510 | we solve it by back-substitution. */ |

511 | |

512 | for (i=p-1; i>=0; i--) { |

513 | j = index[i]; |

514 | z[i] /= Amatrix(i,j); |

515 | for (k=0; k<i; k++) |

516 | z[k] -= Amatrix(k,j) * z[i]; |

517 | } |

518 | |

519 | /* This next part is called the "secondary loop" in the book */ |

520 | while (1) { |

521 | double alpha, tmp; |

522 | |

523 | if (++iter > itmax) /* over the iteration count, oops */ |

524 | goto terminate; |

525 | |

526 | alpha = 2; /* greater than the values we will be |

527 | comparing, which will all be <= 1 */ |

528 | for (i=0; i<p; i++) |

529 | if (z[i] <= 0) { |

530 | tmp = x[index[i]] / (x[index[i]]-z[i]); /* This is <= 1 */ |

531 | if (tmp<alpha) { |

532 | alpha = tmp; |

533 | q = i; /* keep track of the minimum index */ |

534 | } |

535 | } |

536 | |

537 | /* Were all z_j>0 for j in P? If so, set x:=z and break out of |

538 | this inner loop, returning to the main loop. */ |

539 | if (alpha > 1) { |

540 | /* the next four lines of code corresponds to lines 244-247 in |

541 | the Fortran 66 version */ |

542 | for (i=0; i<p; i++) { |

543 | x[index[i]] = z[i]; |

544 | } |

545 | break; |

546 | } |

547 | |

548 | /* Interpolate between old x and new z using 0<alpha<=1 */ |

549 | for (i=0; i<p; i++) |

550 | x[index[i]] += alpha * (z[i] - x[index[i]]); |

551 | |

552 | /* Now we have x[index[q]]=0, so modify A and b and the index |

553 | arrays to move index[q] and any other i with x[i]=0 from set |

554 | P to set Z */ |

555 | |

556 | while (q >= 0) { |

557 | i = index[q]; |

558 | x[i] = 0; /* ensure that rounding errors don't creep in */ |

559 | /* now use Givens rotations to fix up A and b */ |

560 | for (jj=q+1; jj<p; jj++) { |

561 | double c, s; |

562 | ii = index[jj]; |

563 | index[jj-1] = ii; |

564 | GivensCalc(Amatrix(jj-1,ii), Amatrix(jj,ii), |

565 | &c, &s, &Amatrix(jj-1,ii)); |

566 | Amatrix(jj,ii) = 0; |

567 | |

568 | for (l=0; l<n; l++) |

569 | if (l != ii) |

570 | GivensApply(c, s, &Amatrix(jj-1,l), &Amatrix(jj,l)); |

571 | GivensApply(c, s, &b[jj-1], &b[jj]); |

572 | } |

573 | /* We've taken one element out of P */ |

574 | p--; |

575 | index[p] = i; |

576 | |

577 | /* Check that the remaining coefficients in set P are |

578 | feasible, that is, >=0. They should be, because of |

579 | the way alpha was determined. If they are not, |

580 | this is due to rounding error, so we set any |

581 | non-positive elements to zero and move them too |

582 | from set P to set Z */ |

583 | for (q=p-1; q>=0; q--) |

584 | if (x[index[q]] <= 0) break; |

585 | /* now if q>=0, then x[index[q]] <= 0, |

586 | and if q<0, we are done */ |

587 | } /* end of while (q >= 0) loop */ |

588 | |

589 | /* OK, all x[i]=0 have now had their indices moved from P |

590 | to Z. We can now solve the least squares problem |

591 | |Az-b| once again, ready for the next secondary loop. |

592 | So we simply copy b into z and solve the triangular |

593 | system, as before. We could reduce the code size a |

594 | fraction by having this code at the beginning of the |

595 | loop, but that would make the loop exit code (checking |

596 | iter) appear in the middle, which would be a tad |

597 | confusing. */ |

598 | |

599 | for (i=0; i<m; i++) |

600 | z[i] = b[i]; |

601 | for (i=p-1; i>=0; i--) { |

602 | j = index[i]; |

603 | z[i] /= Amatrix(i,j); |

604 | for (k=0; k<i; k++) |

605 | z[k] -= Amatrix(k,j) * z[i]; |

606 | } |

607 | } /* end of secondary loop */ |

608 | } /* end of primary loop */ |

609 | |

610 | terminate: |

611 | sum = 0; |

612 | if (p < m) { |

613 | for (i=p; i<m; i++) |

614 | sum += b[i] * b[i]; |

615 | } else { |

616 | /* this is to get the w (dual) vector correct, in case we ever |

617 | decide to give it as an output to this function. */ |

618 | for (i=0; i<n; i++) |

619 | w[i] = 0; |

620 | } |

621 | |

622 | (*rnorm) = sqrt(sum); |

623 | |

624 | if (preserve) |

625 | free(AA[0]); |

626 | |

627 | return (iter > itmax) ? -3 : 0; |

628 | } |

629 | |

630 | |

631 | // Emacs settings: |

632 | // Local variables: |

633 | // mode: c |

634 | // c-basic-offset: 4 |

635 | // tab-width: 4 |

636 | // End: |

**Note:**See TracBrowser for help on using the repository browser.