Changes in sasmodels/kernel_iq.c [70530778:c57ee9e] in sasmodels
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/kernel_iq.c
r70530778 rc57ee9e 27 27 // parameters in the parameter table. 28 28 // CALL_VOLUME(table) : call the form volume function 29 // CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function 29 30 // CALL_IQ(q, table) : call the Iq function for 1D calcs. 30 31 // CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 32 // CALL_FQ(q, F1, F2, table) : call the Fq function for 1D calcs. 33 // CALL_FQ_A(q, F1, F2, table) : call the Iq function with |q| for 2D data. 31 34 // CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 32 35 // CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes … … 36 39 // PROJECTION : equirectangular=1, sinusoidal=2 37 40 // see explore/jitter.py for definitions. 41 38 42 39 43 #ifndef _PAR_BLOCK_ // protected block so we can include this code twice. … … 188 192 QACRotation *rotation, 189 193 double qx, double qy, 190 double *qa b_out, double *qc_out)194 double *qa_out, double *qc_out) 191 195 { 196 const double dqc = rotation->R31*qx + rotation->R32*qy; 192 197 // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 193 const double dqc = rotation->R31*qx + rotation->R32*qy; 194 const double dqab_sq = -dqc*dqc + qx*qx + qy*qy; 195 //*qab_out = sqrt(fabs(dqab_sq)); 196 *qab_out = dqab_sq > 0.0 ? sqrt(dqab_sq) : 0.0; 198 const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 199 200 *qa_out = dqa; 197 201 *qc_out = dqc; 198 202 } … … 270 274 #endif // _QABC_SECTION 271 275 272 273 276 // ==================== KERNEL CODE ======================== 274 277 #define COMPUTE_F1_F2 defined(CALL_FQ) 275 278 kernel 276 279 void KERNEL_NAME( … … 282 285 global const double *q, // nq q values, with padding to boundary 283 286 global double *result, // nq+1 return values, again with padding 284 const double cutoff // cutoff in the dispersity weight product 287 const double cutoff, // cutoff in the dispersity weight product 288 int32_t effective_radius_type // which effective radius to compute 285 289 ) 286 290 { … … 341 345 // version must loop over all q. 342 346 #ifdef USE_OPENCL 343 double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 344 double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 347 #if COMPUTE_F1_F2 348 double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]); 349 double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]); 350 double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]); 351 double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]); 352 double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]); 353 #else 354 double weight_norm = (pd_start == 0 ? 0.0 : result[nq]); 355 double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]); 356 double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]); 357 double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 358 #endif 345 359 #else // !USE_OPENCL 346 double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 360 #if COMPUTE_F1_F2 361 double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]); 362 double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]); 363 double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]); 364 #else 365 double weight_norm = (pd_start == 0 ? 0.0 : result[nq]); 366 double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]); 367 double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]); 368 #endif 347 369 if (pd_start == 0) { 348 370 #ifdef USE_OPENMP 349 371 #pragma omp parallel for 350 372 #endif 351 for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 373 #if COMPUTE_F1_F2 374 // 2*nq for F^2,F pairs 375 for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0; 376 #else 377 for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 378 #endif 352 379 } 353 //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);380 //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_volume, result[0]); 354 381 #endif // !USE_OPENCL 355 382 … … 442 469 // inner loop and defines the macros that use them. 443 470 444 #if defined(CALL_IQ) 445 // unoriented 1D 471 472 #if defined(CALL_FQ) // COMPUTE_F1_F2 is true 473 // unoriented 1D returning <F> and <F^2> 474 double qk; 475 double F1, F2; 476 #define FETCH_Q() do { qk = q[q_index]; } while (0) 477 #define BUILD_ROTATION() do {} while(0) 478 #define APPLY_ROTATION() do {} while(0) 479 #define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table) 480 481 #elif defined(CALL_FQ_A) 482 // unoriented 2D return <F> and <F^2> 483 double qx, qy; 484 double F1, F2; 485 #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 486 #define BUILD_ROTATION() do {} while(0) 487 #define APPLY_ROTATION() do {} while(0) 488 #define CALL_KERNEL() CALL_FQ_A(sqrt(qx*qx+qy*qy),F1,F2,local_values.table) 489 490 #elif defined(CALL_IQ) 491 // unoriented 1D return <F^2> 446 492 double qk; 447 493 #define FETCH_Q() do { qk = q[q_index]; } while (0) 448 494 #define BUILD_ROTATION() do {} while(0) 449 495 #define APPLY_ROTATION() do {} while(0) 450 #define CALL_KERNEL() CALL_IQ(qk, 496 #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 451 497 452 498 #elif defined(CALL_IQ_A) … … 482 528 #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 483 529 #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 530 484 531 #elif defined(CALL_IQ_XY) 485 532 // direct call to qx,qy calculator … … 495 542 // logic in the IQ_AC and IQ_ABC branches. This will become more important 496 543 // if we implement more projections, or more complicated projections. 497 #if defined(CALL_IQ) || defined(CALL_IQ_A) // no orientation 544 #if defined(CALL_IQ) || defined(CALL_IQ_A) || defined(CALL_FQ) || defined(CALL_FQ_A) 545 // no orientation 498 546 #define APPLY_PROJECTION() const double weight=weight0 499 547 #elif defined(CALL_IQ_XY) // pass orientation to the model … … 645 693 // Note: weight==0 must always be excluded 646 694 if (weight > cutoff) { 647 pd_norm += weight * CALL_VOLUME(local_values.table); 695 weighted_volume += weight * CALL_VOLUME(local_values.table); 696 #if COMPUTE_F1_F2 697 weight_norm += weight; 698 #endif 699 if (effective_radius_type != 0) { 700 weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); 701 } 648 702 BUILD_ROTATION(); 649 703 … … 693 747 } 694 748 #else // !MAGNETIC 695 const double scattering = CALL_KERNEL(); 749 #if COMPUTE_F1_F2 750 CALL_KERNEL(); // sets F1 and F2 by reference 751 #else 752 const double scattering = CALL_KERNEL(); 753 #endif 696 754 #endif // !MAGNETIC 697 755 //printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 698 756 699 757 #ifdef USE_OPENCL 700 this_result += weight * scattering; 758 #if COMPUTE_F1_F2 759 this_F2 += weight * F2; 760 this_F1 += weight * F1; 761 #else 762 this_result += weight * scattering; 763 #endif 701 764 #else // !USE_OPENCL 702 result[q_index] += weight * scattering; 765 #if COMPUTE_F1_F2 766 result[2*q_index+0] += weight * F2; 767 result[2*q_index+1] += weight * F1; 768 #else 769 result[q_index] += weight * scattering; 770 #endif 703 771 #endif // !USE_OPENCL 704 772 } 705 773 } 706 774 } 707 708 775 // close nested loops 709 776 ++step; … … 726 793 // Remember the current result and the updated norm. 727 794 #ifdef USE_OPENCL 728 result[q_index] = this_result; 729 if (q_index == 0) result[nq] = pd_norm; 730 //if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 795 #if COMPUTE_F1_F2 796 result[2*q_index+0] = this_F2; 797 result[2*q_index+1] = this_F1; 798 if (q_index == 0) { 799 result[2*nq+0] = weight_norm; 800 result[2*nq+1] = weighted_volume; 801 result[2*nq+2] = weighted_radius; 802 } 803 #else 804 result[q_index] = this_result; 805 if (q_index == 0) { 806 result[nq+0] = weight_norm; 807 result[nq+1] = weighted_volume; 808 result[nq+2] = weighted_radius; 809 } 810 #endif 811 812 //if (q_index == 0) printf("res: %g/%g\n", result[0], weigthed_volume); 731 813 #else // !USE_OPENCL 732 result[nq] = pd_norm; 733 //printf("res: %g/%g\n", result[0], pd_norm); 814 #if COMPUTE_F1_F2 815 result[2*nq] = weight_norm; 816 result[2*nq+1] = weighted_volume; 817 result[2*nq+2] = weighted_radius; 818 #else 819 result[nq] = weight_norm; 820 result[nq+1] = weighted_volume; 821 result[nq+2] = weighted_radius; 822 #endif 823 //printf("res: %g/%g\n", result[0], weighted_volume); 734 824 #endif // !USE_OPENCL 735 825 736 826 // ** clear the macros in preparation for the next kernel ** 827 #undef COMPUTE_F1_F2 737 828 #undef PD_INIT 738 829 #undef PD_OPEN
Note: See TracChangeset
for help on using the changeset viewer.