Changeset 13124 for trunk/psLib/test/math/tap_psSparse.c
- Timestamp:
- May 1, 2007, 6:20:06 PM (19 years ago)
- File:
-
- 1 edited
-
trunk/psLib/test/math/tap_psSparse.c (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/test/math/tap_psSparse.c
r13123 r13124 1 #include <stdio.h>1 #include <stdio.h> 2 2 #include <string.h> 3 3 #include <pslib.h> … … 6 6 #include "pstap.h" 7 7 8 int main (void)8 int main(void) 9 9 { 10 plan_tests(26); 10 psLogSetFormat("HLNM"); 11 psLogSetLevel(PS_LOG_INFO); 12 plan_tests(40); 13 14 15 // test psSparseAlloc() 16 { 17 psMemId id = psMemGetId(); 18 psSparse *matrix = psSparseAlloc(10, 20); 19 ok(matrix != NULL, "psSparse successfully allocated"); 20 skip_start(matrix == NULL, 12, "skipping tests because psSparseAlloc() returned NULL"); 21 ok(matrix->Aij != NULL, "psSparseAlloc() set ->Aij correctly"); 22 ok(matrix->Aij->n == 0, "psSparseAlloc() set ->Aij->n correctly"); 23 ok(matrix->Si != NULL, "psSparseAlloc() set ->Si correctly"); 24 ok(matrix->Si->n == 0, "psSparseAlloc() set ->Si->n correctly"); 25 ok(matrix->Sj != NULL, "psSparseAlloc() set ->Sj correctly"); 26 ok(matrix->Sj->n == 0, "psSparseAlloc() set ->Sj->n correctly"); 27 ok(matrix->Bfj != NULL, "psSparseAlloc() set ->Bfj correctly"); 28 ok(matrix->Bfj->n == 10, "psSparseAlloc() set ->Bfj->n correctly"); 29 ok(matrix->Qii != NULL, "psSparseAlloc() set ->Qii correctly"); 30 ok(matrix->Qii->n == 10, "psSparseAlloc() set ->Qii->n correctly"); 31 ok(matrix->Nelem == 0, "psSparseAlloc() set ->Nelem correctly"); 32 ok(matrix->Nrows == 10, "psSparseAlloc() set ->Nrows correctly"); 33 skip_end(); 34 psFree(matrix); 35 ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); 36 } 37 11 38 12 39 // test psSparseSolve for a simple normal example matrix 13 40 { 14 41 psMemId id = psMemGetId(); 15 16 //solve a normalized matrix equation with psSparseSolve17 18 42 // the basic equation is Ax = b 19 43 … … 22 46 // Ax = b for x using psSparseSolve. compare with the input values for x. 23 47 24 psSparse *matrix = psSparseAlloc (100, 100);48 psSparse *matrix = psSparseAlloc(100, 100); 25 49 ok(matrix != NULL, "psSparse successfully allocated"); 26 50 skip_start(matrix == NULL, 5, "Skipping tests because psSparseAlloc() failed"); 27 51 28 for (int i = 0; i < 100; i++) {29 psSparseMatrixElement (matrix, i, i, 1.0);52 for(int i = 0; i < 100; i++) { 53 psSparseMatrixElement(matrix, i, i, 1.0); 30 54 if (i + 1 < 100) { 31 psSparseMatrixElement (matrix, i + 1, i, 0.1);55 psSparseMatrixElement(matrix, i + 1, i, 0.1); 32 56 } 33 57 } … … 37 61 psSparseResort(matrix); 38 62 39 psVector *xRef = psVectorAlloc (100, PS_TYPE_F32); 40 for (int i = 0; i < 100; i++) 41 { 63 psVector *xRef = psVectorAlloc(100, PS_TYPE_F32); 64 for (int i = 0; i < 100; i++) { 42 65 xRef->data.F32[i] = 1.0; 43 66 } … … 60 83 float dS = 0; 61 84 float dS2 = 0; 62 for (int i = 0; i < 100; i++) 63 { 85 for (int i = 0; i < 100; i++) { 64 86 float dX = xRef->data.F32[i] - xFit->data.F32[i]; 65 // fprintf (stderr, "%5.3f %5.3f %5.3f %5.3f\n", bVec->data.F32[i], xRef->data.F32[i], xFit->data.F32[i], dX);87 // fprintf(stderr, "%5.3f %5.3f %5.3f %5.3f\n", bVec->data.F32[i], xRef->data.F32[i], xFit->data.F32[i], dX); 66 88 dS2 += PS_SQR(dX); 67 89 dS += dX; 68 90 } 69 91 dS /= 100.0; 70 dS2 = sqrt (dS2/100.0 - dS*dS);71 72 is_float_tol (dS2, 0.0, 1e-4, "scatter: %.20f", dS2);73 74 psFree (matrix);75 psFree (xRef);76 psFree (bVec);77 psFree (xFit);78 79 skip_end(); 80 ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");92 dS2 = sqrt(dS2/100.0 - dS*dS); 93 94 is_float_tol(dS2, 0.0, 1e-4, "scatter: %.20f", dS2); 95 96 psFree(matrix); 97 psFree(xRef); 98 psFree(bVec); 99 psFree(xFit); 100 101 skip_end(); 102 ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); 81 103 } 82 104 … … 84 106 { 85 107 psMemId id = psMemGetId(); 86 87 //solve a non-normalized matrix equation with psSparseSolve88 108 // the basic equation is Ax = b 109 89 110 // create a matrix A with diagonals of 1 and a small number of off diagonal elements. 90 111 // construct a vector x, construct the corresponding vector b by multiplication. solve 91 112 // Ax = b for x using psSparseSolve. compare with the input values for x. 92 113 93 psSparse *matrix = psSparseAlloc (100, 100);114 psSparse *matrix = psSparseAlloc(100, 100); 94 115 ok(matrix != NULL, "psSparse successfully allocated"); 95 116 skip_start(matrix == NULL, 5, "Skipping tests because psSparseAlloc() failed"); 96 117 97 118 for (int i = 0; i < 100; i++) { 98 psSparseMatrixElement (matrix, i, i, 5.0);119 psSparseMatrixElement(matrix, i, i, 5.0); 99 120 if (i + 1 < 100) { 100 psSparseMatrixElement (matrix, i + 1, i, 0.1);121 psSparseMatrixElement(matrix, i + 1, i, 0.1); 101 122 } 102 123 } … … 106 127 psSparseResort(matrix); 107 128 108 psVector *xRef = psVectorAlloc (100, PS_TYPE_F32); 109 for (int i = 0; i < 100; i++) 110 { 129 psVector *xRef = psVectorAlloc(100, PS_TYPE_F32); 130 for (int i = 0; i < 100; i++) { 111 131 xRef->data.F32[i] = 1.0; 112 132 } … … 123 143 constraint.paramDelta = 1e8; 124 144 125 // solve for normalization terms (need include local sky?)145 // solve for normalization terms(need include local sky?) 126 146 psVector *xFit = psSparseSolve(NULL, constraint, matrix, 4); 127 147 … … 129 149 float dS = 0; 130 150 float dS2 = 0; 131 for (int i = 0; i < 100; i++) 132 { 151 for (int i = 0; i < 100; i++) { 133 152 float dX = xRef->data.F32[i] - xFit->data.F32[i]; 134 // fprintf (stderr, "%5.3f %5.3f %5.3f %5.3f\n", bVec->data.F32[i], xRef->data.F32[i], xFit->data.F32[i], dX);153 // fprintf(stderr, "%5.3f %5.3f %5.3f %5.3f\n", bVec->data.F32[i], xRef->data.F32[i], xFit->data.F32[i], dX); 135 154 dS2 += PS_SQR(dX); 136 155 dS += dX; 137 156 } 138 157 dS /= 100.0; 139 dS2 = sqrt (dS2/100.0 - dS*dS);140 141 is_float_tol (dS2, 0.0, 1e-4, "scatter: %.20f", dS2);142 143 psFree (matrix);144 psFree (xRef);145 psFree (bVec);146 psFree (xFit);147 148 skip_end(); 149 ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");158 dS2 = sqrt(dS2/100.0 - dS*dS); 159 160 is_float_tol(dS2, 0.0, 1e-4, "scatter: %.20f", dS2); 161 162 psFree(matrix); 163 psFree(xRef); 164 psFree(bVec); 165 psFree(xFit); 166 167 skip_end(); 168 ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); 150 169 } 151 170 … … 164 183 { 165 184 psMemId id = psMemGetId(); 166 167 //solve a simple, small matrix equation 168 169 // the basic equation (Ax = b) is: 185 // the basic equation(Ax = b) is: 170 186 // |S B'||x| = |f| 171 187 // |B Q ||y| = |g| … … 180 196 181 197 // construct the sparse matrix 182 psSparse *matrix = psSparseAlloc (Nrows, Nrows);198 psSparse *matrix = psSparseAlloc(Nrows, Nrows); 183 199 ok(matrix != NULL, "psSparse successfully allocated"); 184 200 skip_start(matrix == NULL, 5, "Skipping tests because psSparseAlloc() failed"); … … 190 206 191 207 // border region has a width of 1: 192 psSparseBorder *border = psSparseBorderAlloc (matrix, Nborder);208 psSparseBorder *border = psSparseBorderAlloc(matrix, Nborder); 193 209 194 210 // construct the B component: 195 psSparseBorderElementB (border, 0, 0, 0.1);196 psSparseBorderElementB (border, 1, 0, 0.2);211 psSparseBorderElementB(border, 0, 0, 0.1); 212 psSparseBorderElementB(border, 1, 0, 0.2); 197 213 198 214 // construct the T component: 199 psSparseBorderElementT (border, 0, 0, 0.5);215 psSparseBorderElementT(border, 0, 0, 0.5); 200 216 201 217 // construct the X and Y vectors: 202 psVector *xRef = psVectorAlloc (Nrows, PS_TYPE_F32);218 psVector *xRef = psVectorAlloc(Nrows, PS_TYPE_F32); 203 219 xRef->data.F32[0] = 1.0; 204 220 xRef->data.F32[1] = 1.0; 205 psVector *yRef = psVectorAlloc (Nborder, PS_TYPE_F32);221 psVector *yRef = psVectorAlloc(Nborder, PS_TYPE_F32); 206 222 yRef->data.F32[0] = 0.5; 207 223 … … 210 226 211 227 // test the support functions: LowerProduct 212 fVec = psSparseBorderLowerProduct (NULL, border, xRef);213 is_float (fVec->n, 1.0, "f dimen: %d", fVec->n);214 is_float (fVec->data.F32[0], 0.3, "f(0): %f", fVec->data.F32[0]);215 psFree (fVec);228 fVec = psSparseBorderLowerProduct(NULL, border, xRef); 229 is_float(fVec->n, 1.0, "f dimen: %d", fVec->n); 230 is_float(fVec->data.F32[0], 0.3, "f(0): %f", fVec->data.F32[0]); 231 psFree(fVec); 216 232 217 233 // test the support functions: Upper Product 218 fVec = psSparseBorderUpperProduct (NULL, border, yRef);219 is_float (fVec->n, 2.0, "f dimen: %d", fVec->n);220 is_float (fVec->data.F32[0], 0.05, "f(0): %f", fVec->data.F32[0]);221 is_float (fVec->data.F32[1], 0.10, "f(0): %f", fVec->data.F32[1]);222 psFree (fVec);234 fVec = psSparseBorderUpperProduct(NULL, border, yRef); 235 is_float(fVec->n, 2.0, "f dimen: %d", fVec->n); 236 is_float(fVec->data.F32[0], 0.05, "f(0): %f", fVec->data.F32[0]); 237 is_float(fVec->data.F32[1], 0.10, "f(0): %f", fVec->data.F32[1]); 238 psFree(fVec); 223 239 224 240 // test the support functions: Square Product 225 fVec = psSparseBorderSquareProduct (NULL, border, yRef);226 is_float (fVec->n, 1.0, "f dimen: %d", fVec->n);227 is_float (fVec->data.F32[0], 0.25, "f(0): %f", fVec->data.F32[0]);228 psFree (fVec);241 fVec = psSparseBorderSquareProduct(NULL, border, yRef); 242 is_float(fVec->n, 1.0, "f dimen: %d", fVec->n); 243 is_float(fVec->data.F32[0], 0.25, "f(0): %f", fVec->data.F32[0]); 244 psFree(fVec); 229 245 230 246 fVec = NULL; 231 247 psVector *gVec = NULL; 232 psSparseBorderMultiply (&fVec, &gVec, border, xRef, yRef);233 is_float (fVec->data.F32[0], 1.15, "f(0): %f", fVec->data.F32[0]);234 is_float (fVec->data.F32[1], 1.20, "f(1): %f", fVec->data.F32[1]);235 is_float (gVec->data.F32[0], 0.55, "g(0): %f", gVec->data.F32[0]);248 psSparseBorderMultiply(&fVec, &gVec, border, xRef, yRef); 249 is_float(fVec->data.F32[0], 1.15, "f(0): %f", fVec->data.F32[0]); 250 is_float(fVec->data.F32[1], 1.20, "f(1): %f", fVec->data.F32[1]); 251 is_float(gVec->data.F32[0], 0.55, "g(0): %f", gVec->data.F32[0]); 236 252 237 253 // supply the fVec and gVec data to the border … … 253 269 psVector *yFit = NULL; 254 270 psSparseBorderSolve(&xFit, &yFit, constraint, border, 4); 255 is_float_tol (xFit->data.F32[0], 1.0, 1e-4, "f(0): %f", xFit->data.F32[0]);256 is_float_tol (xFit->data.F32[1], 1.0, 1e-4, "f(1): %f", xFit->data.F32[1]);257 is_float_tol (yFit->data.F32[0], 0.5, 1e-4, "g(0): %f", yFit->data.F32[0]);258 259 psFree (xFit);260 psFree (yFit);261 psFree (fVec);262 psFree (gVec);263 psFree (xRef);264 psFree (yRef);265 psFree (border);266 psFree (matrix);267 skip_end(); 268 ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");271 is_float_tol(xFit->data.F32[0], 1.0, 1e-4, "f(0): %f", xFit->data.F32[0]); 272 is_float_tol(xFit->data.F32[1], 1.0, 1e-4, "f(1): %f", xFit->data.F32[1]); 273 is_float_tol(yFit->data.F32[0], 0.5, 1e-4, "g(0): %f", yFit->data.F32[0]); 274 275 psFree(xFit); 276 psFree(yFit); 277 psFree(fVec); 278 psFree(gVec); 279 psFree(xRef); 280 psFree(yRef); 281 psFree(border); 282 psFree(matrix); 283 skip_end(); 284 ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); 269 285 } 270 286 … … 272 288 { 273 289 psMemId id = psMemGetId(); 274 275 // solve a simple, small matrix equation 276 277 // the basic equation (Ax = b) is: 290 // the basic equation(Ax = b) is: 278 291 // |S B'||x| = |f| 279 292 // |B Q ||y| = |g| … … 288 301 289 302 // construct the sparse matrix 290 psSparse *matrix = psSparseAlloc (Nrows, Nrows);303 psSparse *matrix = psSparseAlloc(Nrows, Nrows); 291 304 ok(matrix != NULL, "psSparse successfully allocated"); 292 305 skip_start(matrix == NULL, 5, "Skipping tests because psSparseAlloc() failed"); … … 301 314 302 315 // border region has a width of 1: 303 psSparseBorder *border = psSparseBorderAlloc (matrix, Nborder);316 psSparseBorder *border = psSparseBorderAlloc(matrix, Nborder); 304 317 305 318 // construct the B component: 306 for (int i = 0; i < Nrows; i++) 307 { 319 for (int i = 0; i < Nrows; i++) { 308 320 for (int j = 0; j < Nborder; j++) { 309 psSparseBorderElementB (border, i, j, 0.1);321 psSparseBorderElementB(border, i, j, 0.1); 310 322 } 311 323 } … … 313 325 314 326 // construct the Q component: 315 for (int i = 0; i < Nborder; i++) 316 { 327 for (int i = 0; i < Nborder; i++) { 317 328 for (int j = 0; j < Nborder; j++) { 318 psSparseBorderElementT (border, i, j, i+j+2);329 psSparseBorderElementT(border, i, j, i+j+2); 319 330 } 320 331 } 321 332 322 333 // construct the X and Y vectors: 323 psVector *xRef = psVectorAlloc (Nrows, PS_TYPE_F32); 324 for (int i = 0; i < Nrows; i++) 325 { 334 psVector *xRef = psVectorAlloc(Nrows, PS_TYPE_F32); 335 for (int i = 0; i < Nrows; i++) { 326 336 xRef->data.F32[i] = 1.0; 327 337 } 328 psVector *yRef = psVectorAlloc (Nborder, PS_TYPE_F32); 329 for (int i = 0; i < Nborder; i++) 330 { 338 psVector *yRef = psVectorAlloc(Nborder, PS_TYPE_F32); 339 for (int i = 0; i < Nborder; i++) { 331 340 yRef->data.F32[i] = 1.0; 332 341 } … … 335 344 psVector *fVec = NULL; 336 345 psVector *gVec = NULL; 337 psSparseBorderMultiply (&fVec, &gVec, border, xRef, yRef);346 psSparseBorderMultiply(&fVec, &gVec, border, xRef, yRef); 338 347 339 348 // supply the fVec and gVec data to the border … … 351 360 constraint.paramDelta = 1e8; 352 361 353 // solve for normalization terms (need include local sky?)362 // solve for normalization terms(need include local sky?) 354 363 psVector *xFit = NULL; 355 364 psVector *yFit = NULL; … … 359 368 float dS = 0; 360 369 float dS2 = 0; 361 for (int i = 0; i < Nrows; i++) 362 { 370 for (int i = 0; i < Nrows; i++) { 363 371 float dX = xRef->data.F32[i] - xFit->data.F32[i]; 364 // fprintf (stderr, "%5.3f %5.3f %5.3f %5.3f\n", fVec->data.F32[i], xRef->data.F32[i], xFit->data.F32[i], dX);372 // fprintf(stderr, "%5.3f %5.3f %5.3f %5.3f\n", fVec->data.F32[i], xRef->data.F32[i], xFit->data.F32[i], dX); 365 373 dS2 += PS_SQR(dX); 366 374 dS += dX; 367 375 } 368 376 dS /= Nrows; 369 dS2 = sqrt (dS2/Nrows - dS*dS);370 is_float_tol (dS2, 0.0, 1e-4, "scatter: %.20f", dS2);377 dS2 = sqrt(dS2/Nrows - dS*dS); 378 is_float_tol(dS2, 0.0, 1e-4, "scatter: %.20f", dS2); 371 379 372 380 // measure stdev between yFit and yRef 373 381 float dY = yRef->data.F32[0] - yFit->data.F32[0]; 374 // fprintf (stderr, "%5.3f %5.3f %5.3f %5.3f\n", gVec->data.F32[0], yRef->data.F32[0], yFit->data.F32[0], dS); 375 is_float_tol (dY, 0.0, 2e-4, "scatter: %.20f", dY); 376 377 psFree (xRef); 378 psFree (yRef); 379 psFree (xFit); 380 psFree (yFit); 381 psFree (fVec); 382 psFree (gVec); 383 psFree (matrix); 384 psFree (border); 385 386 skip_end(); 387 ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); 388 } 389 390 return exit_status(); 382 // fprintf(stderr, "%5.3f %5.3f %5.3f %5.3f\n", gVec->data.F32[0], yRef->data.F32[0], yFit->data.F32[0], dS); 383 is_float_tol(dY, 0.0, 2e-4, "scatter: %.20f", dY); 384 385 psFree(xRef); 386 psFree(yRef); 387 psFree(xFit); 388 psFree(yFit); 389 psFree(fVec); 390 psFree(gVec); 391 psFree(matrix); 392 psFree(border); 393 394 skip_end(); 395 ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); 396 } 391 397 }
Note:
See TracChangeset
for help on using the changeset viewer.
