Changeset be43e39 in sasmodels for sasmodels/kernel_iq.c
- Timestamp:
- Oct 19, 2018 5:49:12 PM (6 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 31fc4ad
- Parents:
- e44432d (diff), 353e899 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/kernel_iq.c
r70530778 re44432d 26 26 // MAGNETIC_PARS : a comma-separated list of indices to the sld 27 27 // parameters in the parameter table. 28 // CALL_VOLUME(table) : call the form volume function 28 // CALL_VOLUME(form, shell, table) : assign form and shell values 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. … … 270 274 #endif // _QABC_SECTION 271 275 272 273 276 // ==================== KERNEL CODE ======================== 274 275 277 kernel 276 278 void KERNEL_NAME( … … 282 284 global const double *q, // nq q values, with padding to boundary 283 285 global double *result, // nq+1 return values, again with padding 284 const double cutoff // cutoff in the dispersity weight product 286 const double cutoff, // cutoff in the dispersity weight product 287 int32_t effective_radius_type // which effective radius to compute 285 288 ) 286 289 { … … 341 344 // version must loop over all q. 342 345 #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]); 346 #if defined(CALL_FQ) 347 double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]); 348 double weighted_form = (pd_start == 0 ? 0.0 : result[2*nq+1]); 349 double weighted_shell = (pd_start == 0 ? 0.0 : result[2*nq+2]); 350 double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+3]); 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_form = (pd_start == 0 ? 0.0 : result[nq+1]); 356 double weighted_shell = (pd_start == 0 ? 0.0 : result[nq+2]); 357 double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+3]); 358 double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 359 #endif 345 360 #else // !USE_OPENCL 346 double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 361 #if defined(CALL_FQ) 362 double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]); 363 double weighted_form = (pd_start == 0 ? 0.0 : result[2*nq+1]); 364 double weighted_shell = (pd_start == 0 ? 0.0 : result[2*nq+2]); 365 double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+3]); 366 #else 367 double weight_norm = (pd_start == 0 ? 0.0 : result[nq]); 368 double weighted_form = (pd_start == 0 ? 0.0 : result[nq+1]); 369 double weighted_shell = (pd_start == 0 ? 0.0 : result[nq+2]); 370 double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+3]); 371 #endif 347 372 if (pd_start == 0) { 348 373 #ifdef USE_OPENMP 349 374 #pragma omp parallel for 350 375 #endif 351 for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 376 #if defined(CALL_FQ) 377 // 2*nq for F^2,F pairs 378 for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0; 379 #else 380 for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 381 #endif 352 382 } 353 //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);383 //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_shell, result[0]); 354 384 #endif // !USE_OPENCL 355 385 … … 442 472 // inner loop and defines the macros that use them. 443 473 444 #if defined(CALL_IQ) 445 // unoriented 1D 474 475 #if defined(CALL_FQ) // COMPUTE_F1_F2 is true 476 // unoriented 1D returning <F> and <F^2> 477 double qk; 478 double F1, F2; 479 #define FETCH_Q() do { qk = q[q_index]; } while (0) 480 #define BUILD_ROTATION() do {} while(0) 481 #define APPLY_ROTATION() do {} while(0) 482 #define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table) 483 484 #elif defined(CALL_FQ_A) 485 // unoriented 2D return <F> and <F^2> 486 double qx, qy; 487 double F1, F2; 488 #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 489 #define BUILD_ROTATION() do {} while(0) 490 #define APPLY_ROTATION() do {} while(0) 491 #define CALL_KERNEL() CALL_FQ_A(sqrt(qx*qx+qy*qy),F1,F2,local_values.table) 492 493 #elif defined(CALL_IQ) 494 // unoriented 1D return <F^2> 446 495 double qk; 447 496 #define FETCH_Q() do { qk = q[q_index]; } while (0) 448 497 #define BUILD_ROTATION() do {} while(0) 449 498 #define APPLY_ROTATION() do {} while(0) 450 #define CALL_KERNEL() CALL_IQ(qk, 499 #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 451 500 452 501 #elif defined(CALL_IQ_A) … … 482 531 #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 483 532 #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 533 484 534 #elif defined(CALL_IQ_XY) 485 535 // direct call to qx,qy calculator … … 495 545 // logic in the IQ_AC and IQ_ABC branches. This will become more important 496 546 // if we implement more projections, or more complicated projections. 497 #if defined(CALL_IQ) || defined(CALL_IQ_A) // no orientation 547 #if defined(CALL_IQ) || defined(CALL_IQ_A) || defined(CALL_FQ) || defined(CALL_FQ_A) 548 // no orientation 498 549 #define APPLY_PROJECTION() const double weight=weight0 499 550 #elif defined(CALL_IQ_XY) // pass orientation to the model … … 645 696 // Note: weight==0 must always be excluded 646 697 if (weight > cutoff) { 647 pd_norm += weight * CALL_VOLUME(local_values.table); 698 double form, shell; 699 CALL_VOLUME(form, shell, local_values.table); 700 weight_norm += weight; 701 weighted_form += weight * form; 702 weighted_shell += weight * shell; 703 if (effective_radius_type != 0) { 704 weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); 705 } 648 706 BUILD_ROTATION(); 649 707 … … 693 751 } 694 752 #else // !MAGNETIC 695 const double scattering = CALL_KERNEL(); 753 #if defined(CALL_FQ) 754 CALL_KERNEL(); // sets F1 and F2 by reference 755 #else 756 const double scattering = CALL_KERNEL(); 757 #endif 696 758 #endif // !MAGNETIC 697 759 //printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 698 760 699 761 #ifdef USE_OPENCL 700 this_result += weight * scattering; 762 #if defined(CALL_FQ) 763 this_F2 += weight * F2; 764 this_F1 += weight * F1; 765 #else 766 this_result += weight * scattering; 767 #endif 701 768 #else // !USE_OPENCL 702 result[q_index] += weight * scattering; 769 #if defined(CALL_FQ) 770 result[2*q_index+0] += weight * F2; 771 result[2*q_index+1] += weight * F1; 772 #else 773 result[q_index] += weight * scattering; 774 #endif 703 775 #endif // !USE_OPENCL 704 776 } 705 777 } 706 778 } 707 708 779 // close nested loops 709 780 ++step; … … 726 797 // Remember the current result and the updated norm. 727 798 #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); 799 #if defined(CALL_FQ) 800 result[2*q_index+0] = this_F2; 801 result[2*q_index+1] = this_F1; 802 if (q_index == 0) { 803 result[2*nq+0] = weight_norm; 804 result[2*nq+1] = weighted_form; 805 result[2*nq+3] = weighted_shell; 806 result[2*nq+3] = weighted_radius; 807 } 808 #else 809 result[q_index] = this_result; 810 if (q_index == 0) { 811 result[nq+0] = weight_norm; 812 result[nq+1] = weighted_form; 813 result[nq+2] = weighted_shell; 814 result[nq+3] = weighted_radius; 815 } 816 #endif 817 818 //if (q_index == 0) printf("res: %g/%g\n", result[0], weighted_shell); 731 819 #else // !USE_OPENCL 732 result[nq] = pd_norm; 733 //printf("res: %g/%g\n", result[0], pd_norm); 820 #if defined(CALL_FQ) 821 result[2*nq] = weight_norm; 822 result[2*nq+1] = weighted_form; 823 result[2*nq+2] = weighted_shell; 824 result[2*nq+3] = weighted_radius; 825 #else 826 result[nq] = weight_norm; 827 result[nq+1] = weighted_form; 828 result[nq+2] = weighted_shell; 829 result[nq+3] = weighted_radius; 830 #endif 831 //printf("res: %g/%g\n", result[0], weighted_shell); 734 832 #endif // !USE_OPENCL 735 833
Note: See TracChangeset
for help on using the changeset viewer.