[37] | 1 | ///////////////////////////////////////////////////////////////////////////////// |
---|
| 2 | //// |
---|
| 3 | //// Visualization of the nonzero pattern of JtJ in EPS format. |
---|
| 4 | //// This is part of the sba package, |
---|
| 5 | //// Copyright (C) 2004 Manolis Lourakis (lourakis@ics.forth.gr) |
---|
| 6 | //// Institute of Computer Science, Foundation for Research & Technology - Hellas |
---|
| 7 | //// Heraklion, Crete, Greece. |
---|
| 8 | //// |
---|
| 9 | //// This program is free software; you can redistribute it and/or modify |
---|
| 10 | //// it under the terms of the GNU General Public License as published by |
---|
| 11 | //// the Free Software Foundation; either version 2 of the License, or |
---|
| 12 | //// (at your option) any later version. |
---|
| 13 | //// |
---|
| 14 | //// This program is distributed in the hope that it will be useful, |
---|
| 15 | //// but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 16 | //// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
| 17 | //// GNU General Public License for more details. |
---|
| 18 | //// |
---|
| 19 | /////////////////////////////////////////////////////////////////////////////////// |
---|
| 20 | |
---|
| 21 | #include <stdio.h> |
---|
| 22 | #include <stdlib.h> |
---|
| 23 | #include <math.h> |
---|
| 24 | #include <time.h> |
---|
| 25 | |
---|
| 26 | #include "sba.h" |
---|
| 27 | |
---|
| 28 | #if 0 |
---|
| 29 | #define HORZ_OFFSET 3 // horizontal offset from left edge of paper in PostScript points (~ 1 mm) |
---|
| 30 | #define VERT_OFFSET 6 // vertical offset from lower edge of paper in PostScript points (~ 2 mm) |
---|
| 31 | #else |
---|
| 32 | #define HORZ_OFFSET 0 |
---|
| 33 | #define VERT_OFFSET 0 |
---|
| 34 | #endif |
---|
| 35 | |
---|
| 36 | #define A4_WIDTH 595 // A4 paper width in PostScript points (210 mm, 1 pt = 25.4/72 mm) |
---|
| 37 | //#define A4_HEIGHT 842 // A4 paper height in PostScript points (297 mm) |
---|
| 38 | |
---|
| 39 | /* 1.1x1.1 looks better than 1.0x1.0 */ |
---|
| 40 | #define DOT_WIDTH 1.1 |
---|
| 41 | #define DOT_HEIGHT 1.1 |
---|
| 42 | |
---|
| 43 | /* saves the nonzero pattern of JtJ in EPS format. If mcon!=0, then all U_j and W_ij with j<mcon |
---|
| 44 | * are assumed to be zero and are therefore clipped off. |
---|
| 45 | */ |
---|
| 46 | void sba_hessian2eps(int n, int m, int mcon, int cnp, int pnp, struct sba_crsm *idxij, int *rcidxs, int *rcsubs, |
---|
| 47 | double *U, double *V, double *W, char *fname) |
---|
| 48 | { |
---|
| 49 | register int i, j, ii, jj; |
---|
| 50 | int i0, j0; |
---|
| 51 | int mmcon, nvars, Usz, Vsz, Wsz, nnz, bbox=1, denseblocks=0; |
---|
| 52 | register double *ptr1, *ptr2; |
---|
| 53 | FILE *fp; |
---|
| 54 | time_t tim; |
---|
| 55 | |
---|
| 56 | mmcon=m-mcon; |
---|
| 57 | nvars=mmcon*cnp + n*pnp; |
---|
| 58 | Usz=cnp * cnp; Vsz=pnp * pnp; Wsz=cnp * pnp; |
---|
| 59 | |
---|
| 60 | if((fp=fopen(fname, "w"))==NULL){ |
---|
| 61 | fprintf(stderr, "JtJ2EPS(): failed to open file %s for writing\n", fname); |
---|
| 62 | exit(1); |
---|
| 63 | } |
---|
| 64 | |
---|
| 65 | /* print EPS preamble */ |
---|
| 66 | time(&tim); |
---|
| 67 | fprintf(fp, "%%!PS-Adobe-2.0 EPSF-2.0\n%%%%Title: %s\n%%%%Creator: sba ver. %s\n", fname, SBA_VERSION); |
---|
| 68 | fprintf(fp, "%%%%CreationDate: %s", ctime(&tim)); |
---|
| 69 | fprintf(fp, "%%%%BoundingBox: 0 0 %d %d\n%%%%Magnification: 1.0000\n%%%%Page: 1 1\n%%%%EndComments\n\n", A4_WIDTH, A4_WIDTH+VERT_OFFSET); |
---|
| 70 | fprintf(fp, "/origstate save def\ngsave\n0 setgray\n"); |
---|
| 71 | if(bbox) fprintf(fp, "0 0 %d %d rectstroke\n", A4_WIDTH, A4_WIDTH+VERT_OFFSET); |
---|
| 72 | |
---|
| 73 | /* move the coordinate system to the upper left corner with axes pointing as shown below: |
---|
| 74 | * +------> y |
---|
| 75 | * | |
---|
| 76 | * | |
---|
| 77 | * | |
---|
| 78 | * v x |
---|
| 79 | */ |
---|
| 80 | fprintf(fp, "%d %d translate\n-90 rotate\n", HORZ_OFFSET, A4_WIDTH+VERT_OFFSET); |
---|
| 81 | /* scale the coordinate system so that the hessian pattern fits in the narrowest page dimension */ |
---|
| 82 | fprintf(fp, "%.4lf %.4lf scale\n", ((double)(A4_WIDTH-HORZ_OFFSET))/nvars, ((double)(A4_WIDTH-HORZ_OFFSET))/nvars); |
---|
| 83 | |
---|
| 84 | /* define dot dimensions */ |
---|
| 85 | fprintf(fp, "/w %g def\n/h %g def\n", DOT_WIDTH, DOT_HEIGHT); |
---|
| 86 | /* define shorthand for rectfill */ |
---|
| 87 | fprintf(fp, "/R { rectfill } bind def\n\n"); |
---|
| 88 | |
---|
| 89 | /* process block row j: [0, ..., 0, U_j, 0, ..., 0, W1j, ..., Wnj] */ |
---|
| 90 | for(j=mcon; j<m; ++j){ |
---|
| 91 | ptr1=U + j*Usz; // set ptr1 to point to U_j |
---|
| 92 | nnz=sba_crsm_col_elmidxs(idxij, j, rcidxs, rcsubs); /* find nonzero W_ij, i=0...n-1 */ |
---|
| 93 | |
---|
| 94 | i0=(j-mcon)*cnp; |
---|
| 95 | for(ii=0; ii<cnp; ++ii){ |
---|
| 96 | j0=(j-mcon)*cnp; |
---|
| 97 | for(jj=0; jj<cnp; ++jj){ // row ii of U_j |
---|
| 98 | if(denseblocks || (!denseblocks && ptr1[ii*cnp+jj]!=0.0)) |
---|
| 99 | fprintf(fp, "%d %d w h R\n", i0+ii, j0+jj); |
---|
| 100 | } |
---|
| 101 | |
---|
| 102 | for(i=0; i<nnz; ++i){ |
---|
| 103 | /* set ptr2 to point to W_ij, actual row number in rcsubs[i] */ |
---|
| 104 | ptr2=W + idxij->val[rcidxs[i]]*Wsz; |
---|
| 105 | j0=mmcon*cnp+rcsubs[i]*pnp; |
---|
| 106 | for(jj=0; jj<pnp; ++jj){ // row ii of W_ij |
---|
| 107 | if(denseblocks || (!denseblocks && ptr2[ii*pnp+jj]!=0.0)) |
---|
| 108 | fprintf(fp, "%d %d w h R\n", i0+ii, j0+jj); |
---|
| 109 | } |
---|
| 110 | } |
---|
| 111 | } |
---|
| 112 | } |
---|
| 113 | |
---|
| 114 | /* process block row mmcon+i: [W_i1^T, ..., W_im^T, 0, ..., 0, V_i, 0, ..., 0] */ |
---|
| 115 | for(i=0; i<n; ++i){ |
---|
| 116 | nnz=sba_crsm_row_elmidxs(idxij, i, rcidxs, rcsubs); /* find nonzero W_ij, j=0...m-1 */ |
---|
| 117 | ptr1=V + i*Vsz; // set ptr1 to point to V_i |
---|
| 118 | |
---|
| 119 | i0=cnp*mmcon+i*pnp; |
---|
| 120 | for(ii=0; ii<pnp; ++ii){ |
---|
| 121 | for(j=0; j<nnz; ++j){ |
---|
| 122 | /* set ptr2 to point to W_ij, actual column number in rcsubs[j] */ |
---|
| 123 | if(rcsubs[j]<mcon) continue; /* W_ij is zero */ |
---|
| 124 | |
---|
| 125 | ptr2=W + idxij->val[rcidxs[j]]*Wsz; |
---|
| 126 | |
---|
| 127 | j0=(rcsubs[j]-mcon)*cnp; |
---|
| 128 | for(jj=0; jj<cnp; ++jj){ // row ii of W_ij^T |
---|
| 129 | if(denseblocks || (!denseblocks && ptr2[jj*pnp+ii]!=0.0)) |
---|
| 130 | fprintf(fp, "%d %d w h R\n", i0+ii, j0+jj); |
---|
| 131 | } |
---|
| 132 | } |
---|
| 133 | |
---|
| 134 | j0=cnp*mmcon+i*pnp; |
---|
| 135 | for(jj=0; jj<pnp; ++jj){ // row ii of V_i |
---|
| 136 | if(denseblocks || (!denseblocks && ptr1[ii*pnp+jj]!=0.0)) |
---|
| 137 | fprintf(fp, "%d %d w h R\n", i0+ii, j0+jj); |
---|
| 138 | } |
---|
| 139 | } |
---|
| 140 | } |
---|
| 141 | |
---|
| 142 | fprintf(fp, "grestore\norigstate restore\n\n%%%%Trailer"); |
---|
| 143 | fclose(fp); |
---|
| 144 | } |
---|