Index: unk/psphot/src/psImageJpeg.c =================================================================== --- /trunk/psphot/src/psImageJpeg.c (revision 6530) +++ (revision ) @@ -1,165 +1,0 @@ -# include -# include // for strcasecmp -# include // for unlink -# include -# include "jpeglib.h" - -/* vectors to hold JPEG colormap */ -static psVector *colorMapR = NULL; -static psVector *colorMapG = NULL; -static psVector *colorMapB = NULL; - -# define RANGELIM(A)(PS_MAX(0,PS_MIN(255,(A)))) - -# define SCALEVALUE(VALUE,ZERO,SCALE)(PS_MAX(0,PS_MIN(255,(SCALE*(VALUE-ZERO))))) - -bool psImageJpegColormap (char *name) { - - colorMapR = psVectorRecycle (colorMapR, 256, PS_TYPE_U8); - colorMapG = psVectorRecycle (colorMapG, 256, PS_TYPE_U8); - colorMapB = psVectorRecycle (colorMapB, 256, PS_TYPE_U8); - - /* grayscale */ - if ((!strcasecmp (name, "grayscale")) || (!strcasecmp (name, "greyscale"))) { - for (int i = 0; i < colorMapR->n; i++) { - colorMapR->data.U8[i] = RANGELIM(i); - colorMapG->data.U8[i] = RANGELIM(i); - colorMapB->data.U8[i] = RANGELIM(i); - } - return true; - } - - /* -grayscale */ - if ((!strcasecmp (name, "-grayscale")) || (!strcasecmp (name, "-greyscale"))) { - for (int i = 0; i < colorMapR->n; i++) { - colorMapR->data.U8[i] = RANGELIM(256 - i); - colorMapG->data.U8[i] = RANGELIM(256 - i); - colorMapB->data.U8[i] = RANGELIM(256 - i); - } - return true; - } - - /* rainbow */ - if (!strcasecmp (name, "rainbow")) { - int I1 = 0.25*colorMapR->n; - int I2 = 0.50*colorMapR->n; - int I3 = 0.75*colorMapR->n; - for (int i = 0; i < I1; i++) { - colorMapR->data.U8[i] = 0; - colorMapG->data.U8[i] = 0; - colorMapB->data.U8[i] = RANGELIM(4*i); - } - for (int i = I1; i < I2; i++) { - colorMapR->data.U8[i] = RANGELIM(4*(i - I1)); - colorMapG->data.U8[i] = 0; - colorMapB->data.U8[i] = RANGELIM(4*(I2 - i)); - } - for (int i = I2; i < I3; i++) { - colorMapR->data.U8[i] = 255; - colorMapG->data.U8[i] = 4*(i - I2); - colorMapB->data.U8[i] = 0; - } - for (int i = I3; i < colorMapR->n; i++) { - colorMapR->data.U8[i] = 255; - colorMapG->data.U8[i] = 255; - colorMapB->data.U8[i] = RANGELIM(4*(i - I3)); - } - return true; - } - - /* heat */ - if (!strcasecmp (name, "heat")) { - int I1 = 0.25*colorMapR->n; - int I2 = 0.50*colorMapR->n; - int I3 = 0.75*colorMapR->n; - for (int i = 0; i < I1; i++) { - colorMapR->data.U8[i] = RANGELIM(2*i); - colorMapG->data.U8[i] = 0; - colorMapB->data.U8[i] = 0; - } - for (int i = I1; i < I2; i++) { - colorMapR->data.U8[i] = RANGELIM(2*i); - colorMapG->data.U8[i] = RANGELIM(2*(i - I1)); - colorMapB->data.U8[i] = 0; - } - for (int i = I2; i < I3; i++) { - colorMapR->data.U8[i] = 255; - colorMapG->data.U8[i] = RANGELIM(2*(i - I1)); - colorMapB->data.U8[i] = RANGELIM(2*(i - I2)); - } - for (int i = I3; i < colorMapR->n; i++) { - colorMapR->data.U8[i] = 255; - colorMapG->data.U8[i] = 255; - colorMapB->data.U8[i] = RANGELIM(2*(i - I2)); - } - return true; - } - return false; -} - -bool psImageJpeg (psImage *image, char *filename, float min, float max) { - - struct jpeg_compress_struct cinfo; - struct jpeg_error_mgr jerr; - - int pixel; - JSAMPLE *jpegLine; // Points to data for current line - JSAMPROW jpegLineList[1]; // pointer to JSAMPLE row[s] - JSAMPLE *outPix; - - if (colorMapR == NULL) return false; - if (colorMapG == NULL) return false; - if (colorMapB == NULL) return false; - - /* JPEG init calls */ - cinfo.err = jpeg_std_error (&jerr); - jpeg_create_compress (&cinfo); - - /* open file, prep for jpeg */ - FILE *f = fopen (filename, "w"); - if (f == NULL) { - fprintf (stderr, "failed to open %s for output\n", filename); - return (TRUE); - } - jpeg_stdio_dest(&cinfo, f); - - /* set up color jpeg buffers */ - int quality = 75; - cinfo.image_width = image->numCols; // image width and height, in pixels - cinfo.image_height = image->numRows; - cinfo.input_components = 3; - cinfo.in_color_space = JCS_RGB; - jpeg_set_defaults (&cinfo); - jpeg_set_quality (&cinfo, quality, true); // limit to baseline-JPEG values - jpeg_start_compress (&cinfo, true); - - jpegLine = psAlloc (3*image->numCols*sizeof(JSAMPLE)); - jpegLineList[0] = jpegLine; - - psU8 *Rpix = colorMapR->data.U8; - psU8 *Gpix = colorMapG->data.U8; - psU8 *Bpix = colorMapB->data.U8; - - float zero = min; - float scale = 256.0/(max - min); - - for (int j = 0; j < image->numRows; j++) { - psF32 *row = image->data.F32[j]; - - outPix = jpegLine; - for (int i = 0; i < image->numCols; i++, outPix += 3) { - pixel = SCALEVALUE(row[i],zero,scale); - outPix[0] = Rpix[pixel]; - outPix[1] = Gpix[pixel]; - outPix[2] = Bpix[pixel]; - } - jpeg_write_scanlines (&cinfo, jpegLineList, 1); - } - - jpeg_finish_compress (&cinfo); - fclose (f); - jpeg_destroy_compress (&cinfo); - - psFree (jpegLine); - return true; -} Index: unk/psphot/src/psLibUtils.h =================================================================== --- /trunk/psphot/src/psLibUtils.h (revision 6530) +++ (revision ) @@ -1,18 +1,0 @@ - -# ifndef PS_LIB_UTILS -# define PS_LIB_UTILS - -// structure to carry a dynamic string -typedef struct { - int NLINE; - int Nline; - char *line; -} psLine; - -// psLine functions -- keep out for now? -psLine *psLineAlloc (int Nline); -bool psLineInit (psLine *line); -bool psLineAdd (psLine *line, char *format, ...); -int psStringStrip (char *string); - -# endif Index: unk/psphot/src/psLine.c =================================================================== --- /trunk/psphot/src/psLine.c (revision 6530) +++ (revision ) @@ -1,66 +1,0 @@ -# include -# include "psLibUtils.h" - -static void psLineFree (psLine *line) { - - if (line == NULL) return; - - psFree (line->line); - return; -} - -// allocate a psLine structrue -psLine *psLineAlloc (int Nline) { - - psLine *line; - line = psAlloc (sizeof(psLine)); - psMemSetDeallocator(line, (psFreeFunc) psLineFree); - - line->Nline = 0; - line->NLINE = Nline; - line->line = psAlloc (Nline); - return (line); -} - -bool psLineInit (psLine *line) { - if (line == NULL) return (false); - line->Nline = 0; - return (true); -} - -bool psLineAdd (psLine *line, char *format, ...) { - - int Nchar; - va_list ap; - - if (line == NULL) return (false); - - int nMax = line->NLINE - line->Nline; - - va_start (ap, format); - Nchar = vsnprintf (&line->line[line->Nline], nMax, format, ap); - line->Nline += PS_MIN (nMax - 1, Nchar); - va_end (ap); - - if (Nchar >= nMax) return (false); - return (true); -} - -#ifndef whitespace -#define whitespace(c) (((c) == ' ') || ((c) == '\t')) -#endif - -/* Strip whitespace from the start and end of STRING. */ -int psStringStrip (char *string) { - - int i; - - if (string == (char *) NULL) return (FALSE); - - for (i = 0; whitespace (string[i]); i++); - if (i) memmove (string, string + i, strlen(string+i)+1); - for (i = strlen (string) - 1; (i > 0) && whitespace (string[i]); i--); - string[++i] = 0; - return (i); - -} Index: unk/psphot/src/psPolynomialUtils.c =================================================================== --- /trunk/psphot/src/psPolynomialUtils.c (revision 6530) +++ (revision ) @@ -1,168 +1,0 @@ -# include -# include // for strcasecmp -# include // for unlink -# include - -psPolynomial4D *psVectorChiClipFitPolynomial4D( - psPolynomial4D *poly, - psStats *stats, - const psVector *mask, - psMaskType maskValue, - const psVector *f, - const psVector *fErr, - const psVector *x, - const psVector *y, - const psVector *z, - const psVector *t) -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(z, NULL); - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(t, NULL); - PS_ASSERT_VECTOR_NON_NULL(fErr, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); - - // clipping range defined by min and max and/or clipSigma - float minClipSigma; - float maxClipSigma; - if (isfinite(stats->max)) { - maxClipSigma = +fabs(stats->max); - } else { - maxClipSigma = +fabs(stats->clipSigma); - } - if (isfinite(stats->min)) { - minClipSigma = -fabs(stats->min); - } else { - minClipSigma = -fabs(stats->clipSigma); - } - psVector *fit = NULL; - psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64); - - // eventual expansion: user supplies one of various stats option pairs, - // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to - // evaluate the clipping sigma - // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - - for (int N = 0; N < stats->clipIter; N++) { - int Nkeep = 0; - - poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t); - fit = psPolynomial4DEvalVector (poly, x, y, z, t); - resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit); - - stats = psVectorStats (stats, resid, NULL, mask, maskValue); - psTrace (".psphot.VectorClipFit", 5, "resid stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); - - // set mask if pts are not valid - // we are masking out any point which is out of range - // recovery is not allowed with this scheme - for (int i = 0; i < resid->n; i++) { - if ((mask != NULL) && (mask->data.U8[i] & maskValue)) { - continue; - } - float sigma = hypot (psVectorGet (fErr, i), stats->sampleStdev); - if (resid->data.F64[i] - stats->sampleMedian > sigma*maxClipSigma) { - if (mask != NULL) { - mask->data.U8[i] |= 0x01; - } - continue; - } - if (resid->data.F64[i] - stats->sampleMedian < sigma*minClipSigma) { - if (mask != NULL) { - mask->data.U8[i] |= 0x01; - } - continue; - } - Nkeep ++; - } - - psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n", - Nkeep, x->n); - - stats->clippedNvalues = Nkeep; - psFree (fit); - } - // Free local temporary variables - psFree (resid); - - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); - return(NULL); - } - return(poly); -} - -psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) { - - int ix = x - image->col0; - int iy = y - image->row0; - - psF32 *Fm = &image->data.F32[iy - 1][ix]; - psF32 *Fo = &image->data.F32[iy + 0][ix]; - psF32 *Fp = &image->data.F32[iy + 1][ix]; - - double Fxm = Fm[-1] + Fo[-1] + Fp[-1]; - double Fxp = Fm[+1] + Fo[+1] + Fp[+1]; - double Fym = Fm[-1] + Fm[+0] + Fm[+1]; - double Fyp = Fp[-1] + Fp[+0] + Fp[+1]; - double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1]; - - psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); - poly->mask[2][2] = 1; - poly->mask[1][2] = 1; - poly->mask[2][1] = 1; - - poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ; - - poly->coeff[1][0] = (Fxp - Fxm)/6.0; - poly->coeff[0][1] = (Fyp - Fym)/6.0; - - poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0; - poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0; - - poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0; - - return (poly); -} - -psPlane psImageBicubeMin (psPolynomial2D *poly) { - - psPlane min; - - min.xErr = min.yErr = 0; - - double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]); - - min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det; - min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det; - return (min); -} - -// may need to extend this to change the keyname in the copy -bool psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key) { - - psMetadataItem *item = psMetadataLookup (in, key); - if (item == NULL) return false; - - psMetadataAddItem (out, item, PS_LIST_TAIL, PS_META_REPLACE); - return true; -} Index: unk/psphot/src/psSparse.c =================================================================== --- /trunk/psphot/src/psSparse.c (revision 6530) +++ (revision ) @@ -1,235 +1,0 @@ -# include "psphot.h" - -void psSparseMatrixTest () { - - // build a sparse matrix - psSparse *sparse = psSparseAlloc (3, 9); - - psSparseMatrixElement (sparse, 0, 0, 3.0); - psSparseMatrixElement (sparse, 1, 1, 2.0); - psSparseMatrixElement (sparse, 2, 2, 1.0); - - psSparseMatrixElement (sparse, 1, 0, 0.1); - psSparseMatrixElement (sparse, 2, 0, -0.1); - - psSparseResort (sparse); - for (int i = 0; i < sparse->Nelem; i++) { - fprintf (stderr, "%d %d %f\n", - sparse->Si->data.S32[i], - sparse->Sj->data.S32[i], - sparse->Aij->data.F32[i]); - } - - psVector *x = psVectorAlloc (3, PS_DATA_F32); - x->data.F32[0] = 3; - x->data.F32[1] = 5; - x->data.F32[2] = 7; - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - - psVector *B = psSparseMatrixTimesVector (NULL, sparse, x); - fprintf (stderr, "B: %f %f %f\n", B->data.F32[0], B->data.F32[1], B->data.F32[2]); - - sparse->Bfj->data.F32[0] = B->data.F32[0]; - sparse->Bfj->data.F32[1] = B->data.F32[1]; - sparse->Bfj->data.F32[2] = B->data.F32[2]; - - psSparseConstraint constraint; - constraint.paramMin = -1e8; - constraint.paramMax = +1e8; - constraint.paramDelta = +1e8; - - x = psSparseSolve (x, constraint, sparse, 0); - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - - x = psSparseSolve (x, constraint, sparse, 1); - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - - x = psSparseSolve (x, constraint, sparse, 2); - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - - x = psSparseSolve (x, constraint, sparse, 3); - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - - x = psSparseSolve (x, constraint, sparse, 4); - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - return; -} - -void psSparseResort (psSparse *sparse) { - - int Nelem = sparse->Nelem; - - psVector *index = psVectorSortIndex (NULL, sparse->Sj); - psVector *Aij = sparse->Aij; - psVector *Si = sparse->Si; - psVector *Sj = sparse->Sj; - - // allocate new temporary vectors - psVector *tAij = psVectorAlloc (Nelem, PS_DATA_F32); - psVector *tSi = psVectorAlloc (Nelem, PS_DATA_S32); - psVector *tSj = psVectorAlloc (Nelem, PS_DATA_S32); - for (int i = 0; i < Nelem; i++) { - int j = index->data.U32[i]; - tAij->data.F32[i] = Aij->data.F32[j]; - tSi->data.S32[i] = Si->data.S32[j]; - tSj->data.S32[i] = Sj->data.S32[j]; - } - psFree (index); - psFree (Aij); - psFree (Si); - psFree (Sj); - - sparse->Aij = tAij; - sparse->Si = tSi; - sparse->Sj = tSj; - return; -} - -// user should only add elements above the diagonal, but we don't check this -void psSparseMatrixElement (psSparse *sparse, int i, int j, float value) { - - int k; - - if (i < j) { - fprintf (stderr, "*** error: subdiagonal element ***\n"); - return; - } - - if (i == j) { - // add to the diagonal - sparse->Qii->data.F32[i] = value; - - // check vectors lengths and extend if needed - if (sparse->Nelem >= sparse->Aij->nalloc) { - psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100); - psVectorRealloc (sparse->Si, sparse->Si->nalloc + 100); - psVectorRealloc (sparse->Sj, sparse->Sj->nalloc + 100); - } - - k = sparse->Nelem; - sparse->Aij->data.F32[k] = value; - sparse->Si->data.S32[k] = i; - sparse->Sj->data.S32[k] = j; - - sparse->Nelem ++; - sparse->Aij->n ++; - sparse->Si->n ++; - sparse->Sj->n ++; - } else { - // check vectors lengths and extend if needed - if (sparse->Nelem >= sparse->Aij->nalloc - 1) { - psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100); - psVectorRealloc (sparse->Si, sparse->Si->nalloc + 100); - psVectorRealloc (sparse->Sj, sparse->Sj->nalloc + 100); - } - - k = sparse->Nelem; - sparse->Aij->data.F32[k] = value; - sparse->Si->data.S32[k] = i; - sparse->Sj->data.S32[k] = j; - k++; - - sparse->Aij->data.F32[k] = value; - sparse->Si->data.S32[k] = j; - sparse->Sj->data.S32[k] = i; - - sparse->Nelem += 2; - sparse->Aij->n += 2; - sparse->Si->n += 2; - sparse->Sj->n += 2; - } - return; -} - -void psSparseVectorElement (psSparse *sparse, int i, float value) { - - sparse->Bfj->data.F32[i] = value; - return; -} - -// multiple A * x -psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector) { - - int i, Nelem; - float F; - - if (output == NULL) { - output = psVectorAlloc (vector->n, PS_DATA_F32); - } - - Nelem = 0; - for (int j = 0; j < vector->n; j++) { - F = 0; - while (matrix->Sj->data.S32[Nelem] == j) { - i = matrix->Si->data.S32[Nelem]; - F += vector->data.F32[i] * matrix->Aij->data.F32[Nelem]; - Nelem++; - } - output->data.F32[j] = F; - } - return (output); -} - -psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter) { - - psF32 dG; - - psVector *Qii = sparse->Qii; - psVector *Bfj = sparse->Bfj; - - guess = psVectorCopy (guess, Bfj, PS_DATA_F32); - - // temporary storage for intermediate results - psVector *dQ = psVectorAlloc (guess->n, PS_DATA_F32); - - for (int j = 0; j < Niter; j++) { - dQ = psSparseMatrixTimesVector (dQ, sparse, guess); - for (int i = 0; i < dQ->n; i++) { - dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i]; - if (fabs (dG) > constraint.paramDelta) { - if (dG > 0) { - dG = +constraint.paramDelta; - } else { - dG = -constraint.paramDelta; - } - } - guess->data.F32[i] -= dG; - guess->data.F32[i] = PS_MAX (guess->data.F32[i], constraint.paramMin); - guess->data.F32[i] = PS_MIN (guess->data.F32[i], constraint.paramMax); - } - } - psFree (dQ); - return (guess); -} - -static void psSparseFree (psSparse *sparse) { - if (sparse == NULL) return; - psFree (sparse->Aij); - psFree (sparse->Bfj); - psFree (sparse->Qii); - psFree (sparse->Si); - psFree (sparse->Sj); - return; -} - -// allocate a sparse matrix container for Nrows, with Nelem slots allocated -psSparse *psSparseAlloc (int Nrows, int Nelem) { - - psSparse *sparse = (psSparse *) psAlloc (sizeof(psSparse)); - sparse->Aij = psVectorAlloc (Nelem, PS_DATA_F32); - sparse->Si = psVectorAlloc (Nelem, PS_DATA_S32); - sparse->Sj = psVectorAlloc (Nelem, PS_DATA_S32); - - sparse->Aij->n = 0; - sparse->Si->n = 0; - sparse->Sj->n = 0; - sparse->Nelem = 0; - - sparse->Bfj = psVectorAlloc (Nrows, PS_DATA_F32); - sparse->Qii = psVectorAlloc (Nrows, PS_DATA_F32); - - sparse->Nrows = Nrows; - psMemSetDeallocator(sparse, (psFreeFunc) psSparseFree); - return (sparse); -} - Index: unk/psphot/src/psSparse.h =================================================================== --- /trunk/psphot/src/psSparse.h (revision 6530) +++ (revision ) @@ -1,24 +1,0 @@ - -typedef struct { - double paramDelta; - double paramMin; - double paramMax; -} psSparseConstraint; - -typedef struct { - psVector *Aij; - psVector *Bfj; - psVector *Qii; - psVector *Si; - psVector *Sj; - int Nelem; - int Nrows; -} psSparse; - -void psSparseMatrixTest (); -void psSparseResort (psSparse *sparse); -void psSparseMatrixElement (psSparse *sparse, int i, int j, float value); -void psSparseVectorElement (psSparse *sparse, int i, float value); -psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector); -psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter); -psSparse *psSparseAlloc (int Nrows, int Nelem);