API tools faq
paste
Advertisement
Guest User

example of gcov

a guest
Sep 10th, 2014
552
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 16.99 KB | None | 0 0
raw download clone embed print report
  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
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!