<commit_before>/*
 *
 */

#include <stdlib.h>
#include <string.h>
#include <assert.h>
//#include <omp.h>
#include "config.h"
#include "vhf/fblas.h"
#define MIN(X,Y)        ((X)<(Y)?(X):(Y))
#define BLK     48
#define BUFBASE 96

typedef struct {
        unsigned int addr;
        unsigned char a;
        unsigned char i;
        char sign;
        char _padding;
} _LinkT;
#define EXTRACT_I(I)    (I.i)
#define EXTRACT_A(I)    (I.a)
#define EXTRACT_SIGN(I) (I.sign)
#define EXTRACT_ADDR(I) (I.addr)

double FCI_t1ci_ms0(double *ci0, double *t1, int bcount,
                    int stra_id, int strb_id,
                    int norb, int na, int nlink, _LinkT *clink_index);

static void compress_link(_LinkT *clink, int *link_index,
                          int norb, int nstr, int nlink)
{
        int i, j, k, a, str1, sign;
        for (k = 0; k < nstr; k++) {
                for (j = 0; j < nlink; j++) {
                        a    = link_index[j*4+0];
                        i    = link_index[j*4+1];
                        str1 = link_index[j*4+2];
                        sign = link_index[j*4+3];
                        clink[j].a = a;
                        clink[j].i = i;
                        clink[j].sign = sign;
                        clink[j].addr = str1;
                }
                clink += nlink;
                link_index += nlink * 4;
        }
}

/*
 * t2[:,i,j,k,l] = E^i_j E^k_l|ci0>
 */
static void rdm4_0b_t2(double *ci0, double *t2,
                       int bcount, int stra_id, int strb_id,
                       int norb, int nstrb, int nlinkb, _LinkT *clink_indexb)
{
        const int nnorb = norb * norb;
        const int n4 = nnorb * nnorb;
        int i, j, k, l, a, sign, str1;
        double *t1 = malloc(sizeof(double) * nstrb * nnorb);
        double *pt1, *pt2;
        _LinkT *tab;

        // form t1 which has beta^+ beta |t1> => target stra_id
        FCI_t1ci_ms0(ci0, t1, nstrb, stra_id, 0,
                     norb, nstrb, nlinkb, clink_indexb);

#pragma omp parallel default(none) \
        shared(t1, t2, bcount, strb_id, norb, nlinkb, clink_indexb), \
        private(i, j, k, l, a, str1, sign, pt1, pt2, tab)
{
#pragma omp for schedule(dynamic, 1) nowait
        for (k = 0; k < bcount; k++) {
                memset(t2+k*n4, 0, sizeof(double)*n4);
                tab = clink_indexb + (strb_id+k) * nlinkb;
                for (j = 0; j < nlinkb; j++) {
                        i    = EXTRACT_I   (tab[j]);
                        a    = EXTRACT_A   (tab[j]);
                        str1 = EXTRACT_ADDR(tab[j]);
                        sign = EXTRACT_SIGN(tab[j]);

                        pt1 = t1 + str1 * nnorb;
                        pt2 = t2 + k * n4 + (i*norb+a)*nnorb;
                        if (sign > 0) {
                                for (l = 0; l < nnorb; l++) {
                                        pt2[l] += pt1[l];
                                }
                        } else {
                                for (l = 0; l < nnorb; l++) {
                                        pt2[l] -= pt1[l];
                                }
                        }
                }
        }
}
        free(t1);
}
/*
 * t2[:,i,j,k,l] = E^i_j E^k_l|ci0>
 */
static void rdm4_a_t2(double *ci0, double *t2,
                      int bcount, int stra_id, int strb_id,
                      int norb, int nstrb, int nlinka, _LinkT *clink_indexa)
{
        const int nnorb = norb * norb;
        const int n4 = nnorb * nnorb;
        int i, j, k, l, a, sign, str1;
        double *t1 = malloc(sizeof(double) * nstrb * nnorb);
        double *pt1, *pt2;
        _LinkT *tab = clink_indexa + stra_id * nlinka;

#pragma omp parallel default(none) \
        shared(ci0, t1, t2, bcount, strb_id, norb, nstrb, nlinka, \
               clink_indexa, tab), \
        private(i, j, k, l, a, str1, sign, pt1, pt2)
{
#pragma omp for schedule(dynamic, 1) nowait
        for (j = 0; j < nlinka; j++) {
                i    = EXTRACT_I   (tab[j]);
                a    = EXTRACT_A   (tab[j]);
                str1 = EXTRACT_ADDR(tab[j]);
                sign = EXTRACT_SIGN(tab[j]);

                // form t1 which has alpha^+ alpha |t1> => target stra_id (through str1)
                FCI_t1ci_ms0(ci0, t1, bcount, str1, strb_id,
                             norb, nstrb, nlinka, clink_indexa);<commit_msg>fix bug in 4pdm, omp shared memory conflicts<commit_after>/*
 *
 */

#include <stdlib.h>
#include <string.h>
#include <assert.h>
//#include <omp.h>
#include "config.h"
#include "vhf/fblas.h"
#define MIN(X,Y)        ((X)<(Y)?(X):(Y))
#define BLK     48
#define BUFBASE 96

typedef struct {
        unsigned int addr;
        unsigned char a;
        unsigned char i;
        char sign;
        char _padding;
} _LinkT;
#define EXTRACT_I(I)    (I.i)
#define EXTRACT_A(I)    (I.a)
#define EXTRACT_SIGN(I) (I.sign)
#define EXTRACT_ADDR(I) (I.addr)

double FCI_t1ci_ms0(double *ci0, double *t1, int bcount,
                    int stra_id, int strb_id,
                    int norb, int na, int nlink, _LinkT *clink_index);

static void compress_link(_LinkT *clink, int *link_index,
                          int norb, int nstr, int nlink)
{
        int i, j, k, a, str1, sign;
        for (k = 0; k < nstr; k++) {
                for (j = 0; j < nlink; j++) {
                        a    = link_index[j*4+0];
                        i    = link_index[j*4+1];
                        str1 = link_index[j*4+2];
                        sign = link_index[j*4+3];
                        clink[j].a = a;
                        clink[j].i = i;
                        clink[j].sign = sign;
                        clink[j].addr = str1;
                }
                clink += nlink;
                link_index += nlink * 4;
        }
}

/*
 * t2[:,i,j,k,l] = E^i_j E^k_l|ci0>
 */
static void rdm4_0b_t2(double *ci0, double *t2,
                       int bcount, int stra_id, int strb_id,
                       int norb, int nstrb, int nlinkb, _LinkT *clink_indexb)
{
        const int nnorb = norb * norb;
        const int n4 = nnorb * nnorb;
        int i, j, k, l, a, sign, str1;
        double *t1 = malloc(sizeof(double) * nstrb * nnorb);
        double *pt1, *pt2;
        _LinkT *tab;

        // form t1 which has beta^+ beta |t1> => target stra_id
        FCI_t1ci_ms0(ci0, t1, nstrb, stra_id, 0,
                     norb, nstrb, nlinkb, clink_indexb);

#pragma omp parallel default(none) \
        shared(t1, t2, bcount, strb_id, norb, nlinkb, clink_indexb), \
        private(i, j, k, l, a, str1, sign, pt1, pt2, tab)
{
#pragma omp for schedule(dynamic, 1) nowait
        for (k = 0; k < bcount; k++) {
                memset(t2+k*n4, 0, sizeof(double)*n4);
                tab = clink_indexb + (strb_id+k) * nlinkb;
                for (j = 0; j < nlinkb; j++) {
                        i    = EXTRACT_I   (tab[j]);
                        a    = EXTRACT_A   (tab[j]);
                        str1 = EXTRACT_ADDR(tab[j]);
                        sign = EXTRACT_SIGN(tab[j]);

                        pt1 = t1 + str1 * nnorb;
                        pt2 = t2 + k * n4 + (i*norb+a)*nnorb;
                        if (sign > 0) {
                                for (l = 0; l < nnorb; l++) {
                                        pt2[l] += pt1[l];
                                }
                        } else {
                                for (l = 0; l < nnorb; l++) {
                                        pt2[l] -= pt1[l];
                                }
                        }
                }
        }
}
        free(t1);
}
/*
 * t2[:,i,j,k,l] = E^i_j E^k_l|ci0>
 */
static void rdm4_a_t2(double *ci0, double *t2,
                      int bcount, int stra_id, int strb_id,
                      int norb, int nstrb, int nlinka, _LinkT *clink_indexa)
{
        const int nnorb = norb * norb;
        const int n4 = nnorb * nnorb;
        int i, j, k, l, a, sign, str1;
        double *t1;
        double *pt1, *pt2;
        _LinkT *tab = clink_indexa + stra_id * nlinka;

#pragma omp parallel default(none) \
        shared(ci0, t2, bcount, strb_id, norb, nstrb, nlinka, \
               clink_indexa, tab), \
        private(i, j, k, l, a, str1, sign, t1, pt1, pt2)
{
        t1 = malloc(sizeof(double) * nstrb * nnorb);
#pragma omp for schedule(dynamic, 1) nowait
        for (j = 0; j < nlinka; j++) {
                i    = EXTRACT_I   (tab[j]);
                a    = EXTRACT_A   (tab[j]);
                str1 = EXTRACT_ADDR(tab[j]);
                sign = EXTRACT_SIGN(tab[j]);

                // form t1 which has alpha^+ alpha |t1> => target stra_id (through str1)
                FCI_t1ci_ms0(ci0, t1, bcount, str1, strb_id,
                             norb, nstrb, nlinka, clink_indexa);
{'fraction_non_alphanumeric': 0.12838834070545538, 'fraction_numerical': 0.02200294884881479, 'mean_word_length': 1.3150433184562877, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7279364', 'n_tokens_mistral': 3274, 'n_tokens_neox': 3148, 'n_words': 797}