Advertisement
Guest User

example of gcov

a guest
Sep 10th, 2014
556
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 16.99 KB | None | 0 0
  1.         -:    0:Source:cq_uni.c
  2.         -:    0:Graph:cq_uni.gcno
  3.         -:    0:Data:cq_uni.gcda
  4.         -:    0:Runs:2
  5.         -:    0:Programs:1
  6.         -:    1:/**************************************************************************
  7.         -:    2: chibigq.c      A simple program to compute the CHI index for an
  8.         -:    3: alpha-quartz crystal lattice.
  9.         -:    4:
  10.         -:    5: mdf    2002-12-04  15:56 UTC   Created
  11.         -:    6: tjl    2003-02-20  Implimented IEEE-80 floating point
  12.         -:    7: changed data type real from double to long double
  13.         -:    8: changed Math Lib function calls to their long
  14.         -:    9: equivalents.
  15.         -:   10: ams    2003-05-21  Modified for alpha-quartz crystal lattice.
  16.         -:   11: tjl    2003-05-21  Eliminated program timer, recompiled.
  17.         -:   12: ams    2004-11-28  Modified for benzil crystal lattice, no hydrogens
  18.         -:   13: amm    2005-11-24  Modified for speed by AM Merritt.
  19.         -:   14: ams    2007-08-02  Modified for full benzil unit cell, 78 atoms
  20.         -:   15: ams    2009-06-20  Modified for P3(2) gamma-glycine, 30 atoms
  21.         -:   16: ams    2009-11-12  Modified for P3(2)21 quartz, 9 atoms
  22.         -:   17: smd    2009-12-19  Rewritten in ANSI C for speed (and ease of ASM'ing later)
  23.         -:   18: smd    2009-12-21  Modified for pthread
  24.         -:   19: ***************************************************************************/
  25.         -:   20:
  26.         -:   21:#include <stdio.h>
  27.         -:   22:#include <stdlib.h>
  28.         -:   23:#include <math.h>
  29.         -:   24://#include <pthread.h>
  30.         -:   25:
  31.         -:   26:#define THREADS 4
  32.         -:   27:
  33.         -:   28:/* 18 digits of precision in a long double */
  34.         -:   29:typedef long double real;
  35.         -:   30:
  36.         -:   31:/* bool type */
  37.         -:   32:typedef char boolean;
  38.         -:   33:#define TRUE 1
  39.         -:   34:#define FALSE 0
  40.         -:   35:
  41.         -:   36:/* jacobi transformation of a symmetric 3x3 matrix */
  42.         -:   37:/* this code ripped from my personal library (making it specific to 3x3
  43.         -:   38:   matrices, with no eigenvectors), which is in turn ripped from
  44.         -:   39:   Numerical Recipies oh so long ago, which was in turn probably ripped
  45.         -:   40:   from EISPACK, whose authors no doubt ripped it Jacobi himself,
  46.         -:   41:   who ultimately extracted it from the intellectual aether in 1846. */
  47.        18:   42:void jacobi3(real **A, real (*dest)[3]) {
  48.         -:   43:    int i, j, its, p, q;
  49.         -:   44:    real sm, g, h, c, s, tau, dp, dq, t, th;
  50.         -:   45:    real d[3];
  51.         -:   46:    real a[3][3];
  52.         -:   47:    real b[3], z[3];
  53.         -:   48:
  54.         -:   49:    /* setup initial result */
  55.        72:   50:    for (i = 0; i < 3; ++i)
  56.       216:   51:        for (j = 0; j < 3; ++j)
  57.       162:   52:            a[i][j] = A[i][j];
  58.         -:   53:
  59.         -:   54:    /* update vectors */
  60.        72:   55:    for (i = 0; i < 3; ++i) {
  61.        54:   56:        b[i] = d[i] = a[i][i];
  62.        54:   57:        z[i] = 0.0;
  63.         -:   58:    }
  64.         -:   59:
  65.         -:   60:    /* sweeps */
  66.        36:   61:    for (its = 1; its < 50; ++its) {
  67.         -:   62:        /* done? */
  68.        36:   63:        sm = 0.0;
  69.       108:   64:        for (i = 0; i < 2; ++i)
  70.       180:   65:            for (j = i + 1; j < 3; ++j)
  71.       108:   66:                sm += fabsl(a[i][j]);
  72.         -:   67:
  73.        36:   68:        if (sm == 0.0) {
  74.        18:   69:            for (i = 0; i < 3; ++i) (*dest)[i] = d[i];
  75.         -:   70:            return;
  76.         -:   71:        }
  77.         -:   72:
  78.         -:   73:        /* nuke the off-diagonals */
  79.        54:   74:        for (p = 0; p < 2; ++p) {
  80.        90:   75:            for (q = p + 1; q < 3; ++q) {
  81.        54:   76:                g = 100.0 * fabsl(a[p][q]);
  82.        54:   77:                dp = d[p];
  83.        54:   78:                dq = d[q];
  84.        54:   79:                h = dq - dp;
  85.        54:   80:                if ((fabsl(h) + g) != fabsl(h)) {
  86.        54:   81:                    th = 0.5 * h / a[p][q];
  87.        54:   82:                    t = 1.0 / (fabsl(th) + sqrtl(1.0 + th * th));
  88.        54:   83:                    if (th < 0.0) t = -t;
  89.     #####:   84:                } else t = a[p][q] / h;
  90.        54:   85:                c = 1.0 / sqrtl(1.0 + t * t);
  91.        54:   86:                s = t * c;
  92.        54:   87:                tau = s / (1.0 + c);
  93.        54:   88:                h = t * a[p][q];
  94.        54:   89:                z[p] -= h;
  95.        54:   90:                z[q] += h;
  96.        54:   91:                d[p] -= h;
  97.        54:   92:                d[q] += h;
  98.        54:   93:                a[p][q] = 0.0;
  99.        54:   94:                for (j = 0; j < p - 1; ++j) {
  100.     #####:   95:                    g = a[j][p];
  101.     #####:   96:                    h = a[j][q];
  102.     #####:   97:                    a[j][p] = g - s * (h + g * tau);
  103.     #####:   98:                    a[j][q] = h + s * (g - h * tau);
  104.         -:   99:                }
  105.         -:  100:
  106.        54:  101:                for (j = p + 1; j < q - 1; ++j) {
  107.     #####:  102:                    g = a[p][j];
  108.     #####:  103:                    h = a[j][q];
  109.     #####:  104:                    a[p][j] = g - s * (h + g * tau);
  110.     #####:  105:                    a[j][q] = h + s * (g - h * tau);
  111.         -:  106:                }
  112.         -:  107:
  113.        72:  108:                for (j = q + 1; j < 3; ++j) {
  114.        18:  109:                    g = a[p][j];
  115.        18:  110:                    h = a[q][j];
  116.        18:  111:                    a[p][j] = g - s * (h + g * tau);
  117.        18:  112:                    a[q][j] = h + s * (g - h * tau);
  118.         -:  113:                }
  119.         -:  114:            }
  120.         -:  115:        }
  121.         -:  116:
  122.        72:  117:        for (p = 0; p < 3; ++p) {
  123.        54:  118:            d[p] = (b[p] += z[p]);
  124.        54:  119:            z[p] = 0.0;
  125.         -:  120:        }
  126.         -:  121:    }
  127.     #####:  122:    printf("oh man!\n");
  128.     #####:  123:    exit(EXIT_FAILURE);
  129.         -:  124:}
  130.         -:  125:
  131.         -:  126:
  132.         -:  127:/* sums of integers from lo to hi (inclusive) raised to powers 0, 1, and 2 */
  133.         -:  128:#define SUM0(lo, hi) ((real) ((hi) - (lo) + 1))
  134.         -:  129:#define SUM1(lo, hi) (0.5 * (real) ((hi) + (lo)) * SUM0((lo), (hi)))
  135.         -:  130:#define _SUM2(l, h, n) (((n) + 1.0) * (l) * (l) + (l) * ((n) + 1.0) * (n) + (n) * ((n) + 1.0) * (2.0 * (n) + 1.0) / 6.0)
  136.         -:  131:#define SUM2(l, h) _SUM2((l), (h), (h) - (l))
  137.         -:  132:
  138.         -:  133:
  139.         -:  134:/* Print usage info */
  140.     #####:  135:void usage() {
  141.     #####:  136:    fprintf(stderr, " Usage: chi base step radius ring \n");
  142.     #####:  137:    fprintf(stderr, "        Where ring is 0-false or 1-true \n");
  143.     #####:  138:}
  144.         -:  139:
  145.         -:  140:
  146.         8:  141:void chiqz(real base, real step, real to, real *ANS_A, real *ANS_B, real *ANS_C) {
  147.         -:  142:    real dd, S0, S1, S2, dz, temp1, temp2, oN, chi;
  148.         -:  143:    real d[3];
  149.         -:  144:    int i, j, k, p, nsp;
  150.         -:  145:    real r;
  151.         -:  146:
  152.         -:  147:    /**************** Lattice *****************/
  153.         -:  148:    int lat_i, lat_lo, lat_hi;
  154.         -:  149:    real lat_X0, lat_X1;
  155.         -:  150:    /* the atom sequencer */
  156.         -:  151:    int lat_ia, lat_ib;
  157.         -:  152:    int lat_ia_lo, lat_ia_hi, lat_ib_lo, lat_ib_hi;
  158.         -:  153:    /* lattice specifics */
  159.         -:  154:    real lat_d0, lat_d1;
  160.         -:  155:    /**************** Lattice *****************/
  161.         -:  156:
  162.         -:  157:    /********** sposes **********/
  163.         -:  158:    real *_N;   /* atom counter (good to 2^63)  long double tjl */
  164.         -:  159:    real **_UX; /* mean */
  165.         -:  160:    real ***_U; /* the mixed tensor */
  166.         -:  161:    /********** sposes **********/
  167.         -:  162:
  168.         -:  163:    /* an alpha-quartz crystal lattice generator */
  169.         -:  164:    /* lat_R == radius of sphere */
  170.         -:  165:    /* lattice specifics */
  171.         8:  166:    const real lat_a = 4.9137;
  172.         8:  167:    const real lat_c = 5.4047;
  173.         8:  168:    const real lat_c0 = -0.50;
  174.         8:  169:    const real lat_c1 = sqrtl(0.750);
  175.         -:  170:
  176.         8:  171:    real lat_atx[9] = {0.4697, 0.0000, -0.4697, 0.4133, -0.2672, -0.1461, 0.2672,  0.1461, -0.4133};
  177.         8:  172:    real lat_aty[9] = {0.0000, 0.4697, -0.4697, 0.2672,  0.1461, -0.4133, 0.4133, -0.2672, -0.1461};
  178.         -:  173:    real lat_atz[9] = {1.0 / 6.0, 5.0 / 6.0, 0.5000, 0.1188 + 1.0 / 6.0, 0.1188 + 5.0 / 6.0, 0.6188,
  179.         8:  174:                        -0.1188 - 1.0 / 6.0, 1.0 / 6.0 - 0.1188, 0.3812};
  180.         -:  175:
  181.         8:  176:    const real lat_R = to; /* radius */
  182.         8:  177:    const real lat_RR = lat_R * lat_R;
  183.         8:  178:    const real lat_aa = lat_a * lat_a;
  184.         -:  179:
  185.         -:  180:    /* fill the spheres */
  186.         8:  181:    nsp = (int) ceill((to - base) / step);
  187.         -:  182:    /* Prepare heap for sposes */
  188.         -:  183:
  189.         8:  184:        _N = (real *) malloc(nsp*sizeof(real));
  190.         8:  185:        _UX = (real **) malloc(nsp*sizeof(real*));
  191.         8:  186:        _U = (real ***) malloc(nsp*sizeof(real**));
  192.        26:  187:        for (i = 0; i < nsp; ++i) {
  193.        18:  188:            _UX[i] = (real *) malloc(3*sizeof(real));
  194.        18:  189:            _U[i] = (real **) malloc(3*sizeof(real*));
  195.        72:  190:            for (j = 0; j < 3; ++j) {
  196.        54:  191:                _U[i][j] = (real *) malloc(3*sizeof(real));
  197.        54:  192:                for (k = 0; k < 3; ++k) _U[i][j][k] = 0;
  198.         -:  193:            }
  199.        18:  194:            _N[i] = 0.0;
  200.         -:  195:        }
  201.         -:  196:        /* rewind the generator */
  202.         8:  197:        lat_i = -1;
  203.         8:  198:        lat_ib = lat_ia = lat_ib_hi = lat_ia_hi = 0;
  204.         -:  199:        for (;;) {
  205.         -:  200:            /* return next line through the lattice */
  206.         -:  201:            /* false is returned if there are no more lines */
  207.  43293500:  202:            if (++lat_ib > lat_ib_hi) {
  208.     67700:  203:                if (++lat_ia > lat_ia_hi) {
  209.        80:  204:                    if (++lat_i >= 9) break; /* end? */
  210.         -:  205:                    /* new range for x */
  211.        72:  206:                    dd = lat_R / (lat_a * lat_c1);
  212.        72:  207:                    lat_ia_hi = (int) floorl(dd - lat_atx[lat_i]);
  213.        72:  208:                    lat_ia = (int) ceill(-dd - lat_atx[lat_i]); /* lat_ia_lo */
  214.         -:  209:                }
  215.         -:  210:                /* precompute coordinate */
  216.         -:  211:                lat_d0 = (real) lat_ia + lat_atx[lat_i],
  217.         -:  212:                /* new range for y */
  218.     67692:  213:                dd = sqrtl(lat_RR / lat_aa - lat_c1 * lat_c1 * lat_d0 * lat_d0);
  219.     67692:  214:                lat_ib_hi = (int) floorl(-lat_d0 * lat_c0 + dd - lat_aty[lat_i]);
  220.     67692:  215:                lat_ib = (int) ceill(-lat_d0 * lat_c0 - dd - lat_aty[lat_i]); /* lat_ib_lo */
  221.         -:  216:            }
  222.         -:  217:            /* precompute coordinate */
  223.  43293492:  218:            lat_d1 = (real) lat_ib + lat_aty[lat_i];
  224.         -:  219:            /* precompute atom x/y position */
  225.  43293492:  220:            lat_X0 = lat_a * lat_d0 + lat_a * lat_d1 * lat_c0;
  226.  43293492:  221:            lat_X1 = lat_a * lat_d1 * lat_c1;
  227.         -:  222:
  228. 140700912:  223:            for (i = nsp - 1; i >= 0; --i) {
  229.  97413888:  224:                if (i + 1 == nsp) {
  230.  43293492:  225:                    r = to;
  231.         -:  226:                } else {
  232.  54120396:  227:                    r = base + step * (real) (i + 1);
  233.         -:  228:                }
  234.         -:  229:                /* given a sphere radius, compute the cell endpoints. return false if none. */
  235.         -:  230:                /* new range for 'z' */
  236.  97413888:  231:                dd = r * r - lat_X0 * lat_X0 - lat_X1 * lat_X1;
  237.  97413888:  232:                if (dd < 0.0) break;
  238.  97407420:  233:                dd = sqrtl(dd) / lat_c;
  239.  97407420:  234:                lat_lo = (int) ceill(-dd - lat_atz[lat_i]);
  240.  97407420:  235:                lat_hi = (int) floorl(dd - lat_atz[lat_i]);
  241.         -:  236:                /* add a line of atoms along the z-axis of the crystal */
  242.         -:  237:                /* fill_Spose */
  243.  97407420:  238:                S0 = SUM0(lat_lo, lat_hi);
  244.  97407420:  239:                S1 = SUM1(lat_lo, lat_hi);
  245.  97407420:  240:                S2 = SUM2(lat_lo, lat_hi);
  246.  97407420:  241:                dz = lat_c * lat_atz[lat_i];
  247.  97407420:  242:                temp1 = S0 * lat_X0;
  248.  97407420:  243:                temp2 = S0 * lat_X1;
  249.  97407420:  244:                _N[i] += S0;
  250.  97407420:  245:                _UX[i][0] += temp1;
  251.  97407420:  246:                _UX[i][1] += temp2;
  252.  97407420:  247:                _UX[i][2] += S0 * dz + S1 * lat_c;
  253.  97407420:  248:                _U[i][0][0] += temp1 * lat_X0;
  254.  97407420:  249:                _U[i][0][1] += temp1 * lat_X1;
  255.  97407420:  250:                _U[i][0][2] += temp1 * dz + S1 * lat_X0 * lat_c;
  256.  97407420:  251:                _U[i][1][1] += temp2 * lat_X1;
  257.  97407420:  252:                _U[i][1][2] += temp2 * dz + S1 * lat_X1 * lat_c;
  258.  97407420:  253:                _U[i][2][2] += S0 * dz * dz + 2.0 * lat_c * dz * S1 + lat_c * lat_c * S2;
  259.         -:  254:            }
  260.         -:  255:        }
  261.         -:  256:        /* emit answers */
  262.        26:  257:        for (i = 0; i < nsp; ++i) {
  263.         -:  258:            /* declare the input done: finalize the U tensor */
  264.         -:  259:            /* centre the U */
  265.        18:  260:            oN = (_N[i] == 0.0) ? 0.0 : 1.0 / _N[i];
  266.        18:  261:            _U[i][0][0] -= oN * _UX[i][0] * _UX[i][0];
  267.        18:  262:            _U[i][1][0] = (_U[i][0][1] -= oN * _UX[i][0] * _UX[i][1]);
  268.        18:  263:            _U[i][2][0] = (_U[i][0][2] -= oN * _UX[i][0] * _UX[i][2]);
  269.        18:  264:            _U[i][1][1] -= oN * _UX[i][1] * _UX[i][1];
  270.        18:  265:            _U[i][2][1] = (_U[i][1][2] -= oN * _UX[i][1] * _UX[i][2]);
  271.        18:  266:            _U[i][2][2] -= oN * _UX[i][2] * _UX[i][2];
  272.         -:  267:            /* compute the CHI for the crystal against its 3D enantiomer */
  273.        18:  268:            p = 0;
  274.        18:  269:            jacobi3(_U[i], &d); /* extract eigenvalues */
  275.        18:  270:            for (j = 1; j < 3; ++j) if (d[j] < d[p]) p = j; /* find the smallest */
  276.        18:  271:            chi = 3.0 * d[p] / (_U[i][0][0] + _U[i][1][1] + _U[i][2][2]); /* the CHI is: */
  277.         -:  272:
  278.         -:  273:            //printf("%10.3Lf %20.14Lg %20.18Lf\n", (i + 1 == nsp) ? to : base + step * (real) (i + 1), _N[i], chi);
  279.        18:  274:            ANS_A[i] = (i + 1 == nsp) ? to : base + step * (real) (i + 1);
  280.        18:  275:            ANS_B[i] = _N[i];
  281.        18:  276:            ANS_C[i] = chi;
  282.         -:  277:        }
  283.         -:  278:
  284.         8:  279:}
  285.         -:  280:
  286.         -:  281:
  287.         2:  282:int main(int argc, char **argv) { /* handle the command line */
  288.         -:  283:    int i;
  289.         -:  284:    real base, to, step;
  290.         -:  285:    boolean ring;
  291.         -:  286:    int idx;
  292.         -:  287:    int nsp;
  293.         -:  288:    /* Units per thread */
  294.         -:  289:    real t_step, t_base;
  295.         -:  290:
  296.         -:  291:    /* Answers */
  297.         -:  292:    real *ANS_A, *ANS_B, *ANS_C;
  298.         -:  293:
  299.         -:  294:    /* enforce the Rules and Regulations */
  300.         2:  295:    if (argc != 5) {
  301.     #####:  296:        usage();
  302.     #####:  297:        return EXIT_FAILURE;
  303.         -:  298:    }
  304.         -:  299:
  305.         2:  300:    base = atof(argv[1]);
  306.         2:  301:    step = atof(argv[2]);
  307.         2:  302:    to = atof(argv[3]);
  308.         2:  303:    ring = (boolean) atoi(argv[4]);
  309.         -:  304:
  310.         2:  305:    printf("\n");
  311.         -:  306:
  312.         -:  307:    /* desired computation */
  313.         2:  308:    if (ring) {
  314.         -:  309:        //chiqz(base, step, to);
  315.         2:  310:        nsp = (int) ceill((to - base) / step);
  316.         2:  311:        ANS_A = (real *) malloc(nsp*sizeof(real));
  317.         2:  312:        ANS_B = (real *) malloc(nsp*sizeof(real));
  318.         2:  313:        ANS_C = (real *) malloc(nsp*sizeof(real));
  319.         2:  314:        t_step = (nsp / THREADS) * step;
  320.         2:  315:        t_base = base;
  321.         -:  316:        /*
  322.         -:  317:        printf("nsp=%d\n", nsp);
  323.         -:  318:        printf("THREADS=%d\n", THREADS);
  324.         -:  319:        printf("step=%10.3Lf\n", step);
  325.         -:  320:        printf("t_step=%10.3Lf\n", t_step);
  326.         -:  321:        */
  327.         -:  322:
  328.         8:  323:        for (i = 0; i < THREADS - 1; i++) {
  329.         -:  324:           /* printf("t#%d: base=%10.3Lf to=%10.3Lf\n", i, t_base, t_base + t_step); */
  330.         6:  325:           idx = (nsp / THREADS)*i;
  331.         6:  326:           chiqz(t_base, step, t_base + t_step, &(ANS_A[idx]), &(ANS_B[idx]), &(ANS_C[idx]));
  332.         6:  327:           t_base += t_step;
  333.         -:  328:        };
  334.         -:  329:        /* printf("t#%d: base=%10.3Lf to=%10.3Lf\n", THREADS, t_base, to); */
  335.         2:  330:        idx = (nsp / THREADS)*i;
  336.         2:  331:        chiqz(t_base, step, to, &(ANS_A[idx]), &(ANS_B[idx]), &(ANS_C[idx]));
  337.         -:  332:
  338.        20:  333:        for (i = 0; i < nsp; i++) {
  339.        18:  334:            printf("%10.3Lf %20.14Lg %20.18Lf\n", ANS_A[i], ANS_B[i], ANS_C[i]);
  340.         -:  335:        }
  341.         -:  336:
  342.         -:  337:    } else { /* just one */
  343.         -:  338:        //chiqz(to - step, step, to);
  344.         -:  339:    }
  345.         2:  340:    return EXIT_SUCCESS;
  346.         -:  341:}
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement